{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook lets you generate the elastic response for a unit disk load for terrestrial hydrology and ice sheet applications. Depending on your area of study, the Greens functions generated here are likely inaccurate (especially if your grid cells contain both ocean/lakes and land): use SPOTL/LOADDEF for better results.\n",
    "\n",
    "Susheel Adusumilli, Scripps Institution of Oceanography\n",
    "\n",
    "Last edit: May 10, 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import scipy.io as sio\n",
    "from scipy.interpolate import interp1d\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the domain of interest:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The radius of the disk is: 13.899365830569842 km\n",
      "The disk is approximated using grid cells of size: 0.05 km\n"
     ]
    }
   ],
   "source": [
    "grid_cell = 0.25 #0.25 degree resolution grid\n",
    "earth_radius = 6371.\n",
    "\n",
    "dx = 0.05\n",
    "r = grid_cell*earth_radius*np.pi/360\n",
    "print('The radius of the disk is: ' + str(r) + ' km')\n",
    "print('The disk is approximated using grid cells of size: ' + str(dx) + ' km')\n",
    "rho_w = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first need the Green's function for a point load. For this, I used LoadDef:\n",
    "https://github.com/hrmartens/LoadDef\n",
    "\n",
    "I ran the code in doc/quick_start.txt but first changed the disk_size value in compute_greens_functions.py to 0.25deg to match the size of the disc load considered here. I also edited the code so that it computed the value of the Green's function at 1e-6 and 1e-5 deg angular distance in addition to the default Farrell spacing (you can change this in compute_greens_functions.py). We want the output for the CM reference frame.\n",
    "\n",
    "I included this output in the data/ folder in the repo.\n",
    "\n",
    "We shall now generate the displacement response for a 1m water load over a disk with the above radius 'r'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We are able to approximate the disk with: 99.98859313763498% accuracy\n"
     ]
    }
   ],
   "source": [
    "cm_PREM_ldef = np.loadtxt('data/cm_PREM.txt',skiprows=5) #Load Green's functions\n",
    "cm_PREM_ldef[0,0] = 1.e-7\n",
    "cm_PREM_ldef[1,0] = 1.e-6\n",
    "km_ld = cm_PREM_ldef[:,0]*earth_radius*np.pi/180;\n",
    "u_m_kg = cm_PREM_ldef[:,1]\n",
    "\n",
    "x_ = np.arange(-r,r,dx)\n",
    "y_ = np.arange(-r,r,dx)\n",
    "\n",
    "x,y = np.meshgrid(x_,y_)\n",
    "mask = np.sqrt(x**2 + y**2) <= r #Approximation of a disk\n",
    "\n",
    "dg_cm = np.arange(0.0001,5001.0001,1) #Go up to 5000 km away from source\n",
    "ug_cm = np.zeros_like(dg_cm)\n",
    "\n",
    "for i in np.arange(0,len(x)):\n",
    "    for j in np.arange(0,len(y)):\n",
    "        if mask[i,j]:\n",
    "            d_=np.sqrt((x[i,j] - dg_cm)**2 + y[i,j]**2)\n",
    "            fi = interp1d(km_ld, u_m_kg, kind ='slinear', fill_value=\"extrapolate\")\n",
    "            ug_cm = ug_cm + rho_w*fi(d_)*(dx*1000)**2\n",
    "\n",
    "approx_factor = (np.pi*r**2)/(np.sum(mask[:])*dx**2)\n",
    "area_disk = (np.sum(mask[:])*dx**2)\n",
    "print('We are able to approximate the disk with: ' + str(100/approx_factor) + '% accuracy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEKCAYAAAAvlUMdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXxU9dX48c+ZmUz2BULCqiyyBghhFdwoUlAEpaJWbRVFraKiba1W+2grLn266E+t1cet4t7SulFc6oKKShERFFlkR1D2hCUQsmfO7487CQGSyc0ySSY579frvmbm3vu9c2YIObnfVVQVY4wxpqF5mjoAY4wxLZMlGGOMMWFhCcYYY0xYWIIxxhgTFpZgjDHGhIUlGGOMMWER1gQjImeKyFoR2SAit1VxXETk4eDx5SIypKayIjJTRLaJyLLgdlY4P4Mxxpi6CVuCEREv8CgwAcgALhaRjKNOmwD0Cm5XA4+5LPugqmYFt7fD9RmMMcbUna+mE0TEAwwCOgEFwCpV3eXi2iOADaq6KXid2cBk4JtK50wGnldntOciEUkRkY5ANxdljTHGNGPVJhgROQG4FfghsB7IBmKA3iKSDzwBPKeqgWou0Rn4vtLrrcCJLs7p7KLsDBGZCiwBfqWq+6qI/2qcuyLi4+OH9u3bt7qPWrX930HRQWjfv3bljDGmhVi6dGmOqqbVtXyoO5h7caqsrtGj5pMRkXTgJ8ClwHPVlJcq9h09L01154Qq+xhwT/D1PcD/A6445mTVJ4EnAYYNG6ZLliypJsxqvP1rWD4bbqtlOWOMaSFEZEt9ylebYFT14hDHdgMP1XDtrcBxlV53Aba7PMdfXdnK1XMi8hTwZg1x1I0/Dorzw3JpY4xpDdy0wXiBiTjtIhXnq+oDNRT9AuglIt2BbcBFOHc9lc3Fqe6ajVMFlquqO0Qku7qyItJRVXcEy58LrKzpM9SJPx4CJVBaDD5/WN7CGGNashoTDPAGUAisAKprbzmGqpaKyAzgXcALzFLVVSIyPXj8ceBt4CxgA5APTAtVNnjpP4tIFk4V2WbgGrcx1UpUvPNYcsgSjDHG1IGbBNNFVTPrcvFgF+K3j9r3eKXnClzvtmxw/6V1iaXW/HHOY3E+xLZplLdsSCUlJWzdupXCwsKmDsUY08zFxMTQpUsXoqKiGvS6bhLMf0RkvKq+16Dv3Nz5E5zHkshsh9m6dSuJiYl069YNkar6TBhjDKgqe/bsYevWrXTv3r1Br+1moOUi4HURKRCRAyJyUEQONGgUzVFU+R1MXtPGUUeFhYWkpqZacjHGhCQipKamhqW2w80dzP8DRgErju6u3KJVriKLUJZcjDFuhOt3hZs7mPXAylaVXKBSI3/kJhhjjGlKbhLMDmC+iPxGRG4q38IdWJPzBxNMhFaRNQder5esrCz69+/PoEGDeOCBBwgEnI6IS5Ys4cYbb6y27Pz585k0aVLI67s5x60f/OAH1HowrjEmJDdVZN8GN39wax1aQBVZU4uNjWXZsmUA7N69m5/85Cfk5uZy1113MWzYMIYNG9bEEdZOWVkZXq+3qcMwJmLUeAejqndVtTVGcE3KqsgaVHp6Ok8++SSPPPIIqnrE3cfHH39MVlYWWVlZDB48mIMHDx5R9osvvmDw4MFs2rSp2uvv3buXH/3oR2RmZjJy5EiWL18OwOLFiznppJMYPHgwJ510EmvXrgWgoKCAiy66iMzMTC688EIKCgqqvG63bt24++67OeWUU3j55Zd57733GDVqFEOGDOGCCy4gL8+5w73tttvIyMggMzOTm2++GYDLL7+c6dOnc+qpp9K7d2/efNOZdKKwsJBp06YxcOBABg8ezEcffQTAs88+y5QpUzjzzDPp1asXv/71rwEnsV1++eUMGDCAgQMH8uCDDwKwceNGzjzzTIYOHcqpp57KmjVrav8PY0wYhZrs8kngr6q6oopj8cCFQJGqvhTG+JpORRXZoaaNowHc9cYqvtnesB3/MjolcefZtZsItEePHgQCAXbv3n3E/vvvv59HH32Uk08+mby8PGJiYiqOLVy4kBtuuIF///vfHH/88dVe+84772Tw4MHMmTOHDz/8kKlTp7Js2TL69u3LJ598gs/nY968efzP//wPr776Ko899hhxcXEsX76c5cuXM2TIkGqvHRMTw4IFC8jJyWHKlCnMmzeP+Ph4/vSnP/HAAw8wY8YMXn/9ddasWYOIsH///oqymzdv5uOPP2bjxo2MGTOGDRs28OijjwKwYsUK1qxZw/jx41m3bh0Ay5Yt46uvviI6Opo+ffpwww03sHv3brZt28bKlc6kFeXXv/rqq3n88cfp1asXn3/+Oddddx0ffvhhrf5NjAmnUFVk/wf8VkQG4kzHUj6bci8gCZgFtMzkAhAVC0iLSDDNSVV9RU4++WRuuukmfvrTnzJlyhS6dOkCwOrVq7n66qt577336NSpU8jrLliwgFdffRWA008/nT179pCbm8uBAwe47LLLWL9+PSJCSUkJAJ988klFG1BmZiaZmdWPJb7wwgsBWLRoEd988w0nn3wyAMXFxYwaNYqkpCRiYmK46qqrmDhx4hHtQj/+8Y/xeDz06tWLHj16sGbNGhYsWMANN9wAQN++fenatWtFghk7dizJyckAZGRksGXLFvr378+mTZu44YYbmDhxIuPHjycvL4+FCxdywQUXVLxXUVFRyO/ImMYWarLLZcCPRSQBGAZ0xFkPZrWqrm2k+JqOiDMWpgVUkdX2TiNcNm3ahNfrJT09ndWrV1fsv+2225g4cSJvv/02I0eOZN68eQB07NiRwsJCvvrqqxoTTFWJS0T47W9/y5gxY3j99dfZvHkzP/jBD4447kZ8fHzFe4wbN45//OMfx5yzePFiPvjgA2bPns0jjzxScSdx9HuISJWxlouOjq547vV6KS0tpU2bNnz99de8++67PProo/zrX//ioYceIiUlpaKNy5jmyE0bTJ6qzlfVf6jqnFaRXMpFJ0BRyx9T2hiys7OZPn06M2bMOOaX7saNGxk4cCC33norw4YNq2hLSElJ4a233uJ//ud/mD9/fsjrn3baabz0knNDPX/+fNq1a0dSUhK5ubl07twZcNo4qjp/5cqVFW02oYwcOZL//ve/bNiwAYD8/HzWrVtHXl4eubm5nHXWWTz00ENH/NJ/+eWXCQQCbNy4kU2bNtGnT58j3nvdunV899139OnTp9r3zcnJIRAIcN5553HPPffw5ZdfkpSURPfu3Xn55ZcBJ/l9/fXXNX4GYxqTm15krVd8GhzKaeooIlZBQQFZWVmUlJTg8/m49NJLuemmY3u4P/TQQ3z00Ud4vV4yMjKYMGECn332GQDt27fnjTfeYMKECcyaNYsTTzx6zTrHzJkzmTZtGpmZmcTFxfHcc84yRb/+9a+57LLLeOCBBzj99NMrzr/22msrzs/KymLEiBE1fp60tDSeffZZLr744orqqHvvvZfExEQmT55MYWEhqlrRCA/Qp08fRo8eza5du3j88ceJiYnhuuuuY/r06QwcOBCfz8ezzz57xJ3L0bZt28a0adMqunj/4Q9/AOCll17i2muv5d5776WkpISLLrqIQYMG1fg5jGks0hrGT9ZpwTGA53/krGr5sw8aPqgwW716Nf369WvqMFq1yy+/nEmTJnH++ec3dSjG1Kiq3xkislRV6zyeIGQVmYh4ReS+ul484iWkw6HdNZ9njDHmGCGryFS1TESGioi0uqliwKrITL1UbvMxpjVy0wbzFfBvEXkZqOizq6qvhS2q5iI+zelFVpTnNPgbY4xxzU2CaQvsAU6vtE+Blp9gEtKdx0O7LcEYY0wt1ZhgVHVaYwTSLMUHE0xeNrTt0bSxGGNMhKlxHIyI9BaRD0RkZfB1pojcEf7QmoGENOfRGvqNMabW3EzX/xTwG6AEQFWXAxeFM6hmo+IOxhJMbe3Zs6diAssOHTrQuXNnsrKySElJISMjo1FjmTNnDt98803F69/97ncVswXUxubNmxkwYEBDhtYiNOSyCVV59tlnmTFjBgCPP/44zz///DHnVPdvs3379lp3E3/ooYfIzz88g0dCQv2qxy+++GIyMzOPGB9Vn2u4+fm9/PLLeeWVV47ZH+5/q6O5aYOJU9XFR42+Lg1TPM1LfDvn8VB208YRgVJTUytGtM+cOZOEhARuvvlmNm/eHJYf8NLSUny+qn+c58yZw6RJkyoS2913393g71/OpvQPr+nTp9fq/E6dOlX5izaUhx56iEsuuYS4uLhalavKzp07WbhwIVu2bHFd5uif5bpco7lwcweTIyIn4DTsIyLn4yxC1vJ5oyC2jSWYBlZWVsbPfvYz+vfvz/jx4yumyq9u+vktW7YwduxYMjMzGTt2LN999x3g/JV20003MWbMGG699dYqyy9cuJC5c+dyyy23kJWVxcaNG4/46+6LL77gpJNOYtCgQYwYMYKDBw+yefNmTj31VIYMGcKQIUNYuHBhyM8zf/58xowZw09+8hMGDhxIWVkZt9xyC8OHDyczM5MnnngCgB07dnDaaaeRlZXFgAED+PTTTwHnL+Rf/epXDBkyhLFjx5Kd7fy8LVu2jJEjR5KZmcm5557Lvn37AGdxtFtvvZURI0bQu3fviuusWrWKESNGkJWVRWZmJuvXrwfgxRdfrNh/zTXXUFZWFvLz1PZ9ywUCAXr16lURfyAQoGfPnuTk5BxxTrdu3Y6Ycbpnz57s2rWLN954gxNPPJHBgwfzwx/+kF27dh0T28yZM7n//vsBWLp0KYMGDWLUqFEVM1QfrfKdTXXLIVT28MMPs337dsaMGcOYMWMq9t9+++0MGjSIkSNHVsSVnZ3Neeedx/Dhwxk+fDj//e9/j7ne+PHj2b17N1lZWXz66ac89dRTDB8+nEGDBnHeeedV3Ckd/bMc6hqVf36XLl3K6NGjGTp0KGeccQY7dhz7q/mdd96hb9++nHLKKbz2WiP3zVLVkBvQA5gH5APbgAVA15rKNadt6NChWmd/Ha46+5K6l28i33zzzeEXb9+qOuusht3evtV1LHfeeafed999qqr67bffqtfr1a+++kpVVS+44AJ94YUXVFX19NNP13Xr1qmq6qJFi3TMmDGqqjpp0iR99tlnVVX16aef1smTJ6uq6mWXXaYTJ07U0tLSkOUvu+wyffnllyviKX9dVFSk3bt318WLF6uqam5urpaUlOihQ4e0oKBAVVXXrVun5T8/3377rfbv3/+Yz/fRRx9pXFycbtq0SVVVn3jiCb3nnntUVbWwsFCHDh2qmzZt0vvvv1/vvfdeVVUtLS3VAwcOqKoqoC+++KKqqt511116/fXXq6rqwIEDdf78+aqq+tvf/lZ//vOfq6rq6NGj9aabblJV1bfeekvHjh2rqqozZsyouE5RUZHm5+frN998o5MmTdLi4mJVVb322mv1ueeeC/nvVdv3/eijj3TixImqqjpz5kx98MEHVVX13Xff1SlTphxz/RtvvFFnzZqlqs6/U/l19u7dq4FAQFVVn3rqqYr3euaZZyq+k8o/S5XjvPnmm6v8t6n8b/bMM89o9+7ddf/+/VpQUKDHH3+8fvfdd8eU6dq1q2ZnZ1e8BnTu3LmqqnrLLbdU/NtefPHF+umnn6qq6pYtW7Rv374h319VNScnp+L57bffrg8//LCqHvuzHOoa5T+/xcXFOmrUKN29e7eqqs6ePVunTZt2xDkFBQXapUsXXbdunQYCAb3gggsq/q2OdsTvjMOffYnW43evmyoyVdUfBteA8ajqQRHpHp501wwlpNsdTAPr3r07WVlZAAwdOpTNmzeHnH7+s88+q/jL69JLLz3iL88LLrgAr9dbp+nr165dS8eOHRk+fDgASUlJABw6dIgZM2awbNkyvF5vxVT6oYwYMYLu3Z3/Fu+99x7Lly+v+CszNzeX9evXM3z4cK644gpKSkr40Y9+VPEdeDyeiiUBLrnkEqZMmUJubi779+9n9OjRAFx22WVHfLYpU6Yc8f0BjBo1it///vds3bqVKVOm0KtXLz744AOWLl1a8RkLCgpIT0+v9nPU5X0ru+KKK5g8eTK/+MUvmDVrFtOmHdsJ9cILL+Tuu+9m2rRpzJ49u+Kzb926lQsvvJAdO3ZQXFxc8X26ifPSSy/lP//5T7Xnl6tqOYTjjjsuZBm/319RrTt06FDef/99AObNm3dE296BAwc4ePAgiYmJ1V5r5cqV3HHHHezfv5+8vDzOOOOMimPlP8turV27lpUrVzJu3DjAqRno2LHjEeesWbOG7t2706tXL8D5+XryySddv0d9uUkwrwJDVLXywiivAEPDE1IzE58GOyJ8ltoJf2zqCI5w9JT0BQUFBAIB19PPV24PLJ9Kvzbly6lqlVP2P/jgg7Rv356vv/6aQCBwxAJo1SmPo/y6f/3rX4/45VHuk08+4a233uLSSy/llltuYerUqcec42YZgfLvsHxKf4Cf/OQnnHjiibz11lucccYZ/O1vf0NVueyyyyomyKyvqt63suOOO4727dvz4Ycf8vnnn1fMGl3ZqFGj2LBhA9nZ2cyZM4c77nA6pd5www3cdNNNnHPOOcyfP5+ZM2dWG0d1/3Zu4w/1GY4WFRVV8V6VywQCAT777DNiY2Ndv//ll1/OnDlzGDRoEM8+++wRs4RX/hlyQ1Xp379/xcSw1anL99RQqm2DEZG+InIekCwiUyptl+MsPNY62B1Mowg1/fxJJ53E7NmzAWcG4VNOOaVW5RMTE49Zhhmcxb62b9/OF198AcDBgwcpLS0lNzeXjh074vF4eOGFF2psszjaGWecwWOPPVaxuNm6des4dOgQW7ZsIT09nZ/97GdceeWVfPnll4Dzi6r8bufvf/87p5xyCsnJybRp06aineOFF16o+Gu9Ops2baJHjx7ceOONnHPOOSxfvpyxY8fyyiuvVKwiunfv3orG4qlTp7J48eIjrlGX9z3aVVddxSWXXMKPf/zjKv8iFxHOPfdcbrrpJvr160dqairAEUsrlM+GXZ2UlBSSk5NZsGABQJWJrK6q+3k52vjx43nkkUcqXrv54+bgwYN07NiRkpKSesfcp08fsrOzKxJMSUkJq1atOuKcvn378u2337Jx40aAKtcyCqdQjfx9gElACnB2pW0I8LPwh9ZMxKc5a8KUFDZ1JC3eSy+9xNNPP82gQYPo378///73vwGn4fWZZ54hMzOTF154gb/85S+1Kn/RRRdx3333MXjw4Ir/aOBUffzzn//khhtuYNCgQYwbN47CwkKuu+46nnvuOUaOHMm6detq/ZflVVddRUZGBkOGDGHAgAFcc801lJaWMn/+fLKyshg8eDCvvvoqP//5zwHnL9dVq1YxdOhQPvzwQ373u98Bzi/ZW265hczMTJYtW1axvzr//Oc/GTBgAFlZWaxZs4apU6eSkZHBvffey/jx48nMzGTcuHEVDcHLly8/pkqlLu97tHPOOYe8vLwqq8fKXXjhhbz44osV1WPgNOBfcMEFnHrqqbRr167G93nmmWe4/vrrGTVqVK3uImpy9dVXM2HChCMa+avy8MMPs2TJEjIzM8nIyODxxx+v8dr33HMPJ554IuPGjaNv3771itPv9/PKK69w6623MmjQILKyso7pkBITE8OTTz7JxIkTOeWUU+jatWu93rO2apyuX0RGqWroe7Bmrs7T9QMsfQ7euBF+sQJSql8Tvrmx6fojR0JCAnl5eY36ngcOHODKK6+suONrSEuWLOGXv/zlMb3MTPMWjun6XU12KSLXA/2pVDWmqlfU9U0jSkKl6WIiKMEYE0pSUlJYkssf//hHHnvssQatsjKRy804mBeADsAZwMdAF6DmCsqWonw0v7XDmDBp7LuXcLrtttvYsmVLle1kpvVxk2B6qupvgUOq+hwwERgY3rCakQiej6ym6k9jjIHw/a5wk2BKgo/7RWQAkAx0C0s0zVGEzkcWExPDnj17LMkYY0JSVfbs2eOqO35tuWmDeVJE2gC/BeYCCcHnrUNUDEQnRVwVWZcuXdi6dWvFtB3GGFOdmJgYunTp0uDXdbMezN+CTz/GmTam9YlPi7g7mKioqJAjoY0xJtxqTDAishFYBHwKfKKq39RQpOWJT4u4OxhjjGlqbtpgMoAngFTgfhHZJCKvhzesZiYh8u5gjDGmqblJMGU4Df1lQADYBbj6bSsiZ4rIWhHZICK3VXFcROTh4PHlIjKkFmVvFhEVkZqH/NZXvE0XY4wxteWmkf8AsAJ4AHhKVfe4ubCIeIFHgXHAVuALEZl7VBXbBKBXcDsReAw4saayInJc8Nh3bmKpt4R0KNgLZSXOGjHGGGNq5OYO5mLgE+A6YLaI3CUiY12UGwFsUNVNqloMzAYmH3XOZOD54NIDi4AUEenoouyDwK8JLoIWdvHlY2FyQp9njDGmQo0JRlX/raq3ANcAbwOXA2+6uHZn4PtKr7cG97k5p9qyInIOsE1VQ86hLyJXi8gSEVlS76665dPFROBgS2OMaSo1JhgReTXYk+wvQDwwFWjj4tpVLUJw9B1HdedUuV9E4oDbgRqnd1XVJ1V1mKoOS0tLqzHYkMrvYPKsHcYYY9xy0wbzR+BLVa3dohjOXUflpeK6ANtdnuOvZv8JQHfg6+AiOl2AL0VkhKrurGV87sVH7nQxxhjTVNxUkX1Rh+QC8AXQS0S6i4gfuAhnJoDK5gJTg73JRgK5qrqjurKqukJV01W1m6p2w0lQQ8KaXKDSjMqWYIwxxi03dzB1oqqlIjIDeBfwArNUdZWITA8efxynTecsYAOQD0wLVTZcsdbInwC+WOuqbIwxtRC2BAOgqm/jJJHK+x6v9FyB692WreKcbvWP0gURZ7ClJRhjjHHNVYIRkUycGZQrzlfV18IUU/MUn25VZMYYUwtu5iKbBWQCq3BG8oPT06t1JZiEdNjfOOM6jTGmJXBzBzNSVTPCHklzF98Oti5p6iiMMSZiuBnJ/5mIWIKJT4f8HAjUpUOdMca0Pm7uYJ7DSTI7gSKcQZCqqplhjay5SUgHDUD+3sPLKBtjjKmWmwQzC7gUZ8LLQA3ntlwVgy2zLcEYY4wLbhLMd6p69ADJ1ueI+cisxtAYY2riJsGsEZG/A2/gVJEBrbSbMth8ZMYY45KbBBOLk1jGV9rXCrsp23xkxhhTGzUmGFWd1hiBNHsxKeCJssGWxhjjkptuyscQkRqny29xRJyGfpsuxhhjXKlTggGuatAoIkVCmt3BGGOMS9VWkYnIgeoO4bTLtD7x6XYHY4wxLoVqg9kPDFfVXUcfEJHvqzi/5UtIh93fNHUUxhgTEUJVkT0PdK3m2N/DEEvzV94Go0ev/GyMMeZo1d7BqOodIY7dGp5wmrmEdCgrhsL9ENumqaMxxphmLawLjjUXq7YfYNi97xMf7SM5NoquqfF0bxdPj3bxZB2XQtfUOESk5guVTxeTl20JxhhjatAqEkzbeD/j+3fgUFEpew8Vs+z7fby1fDuBYE1XemI0I7q3ZXTvNM4c0IHEmKiqLxRfabBlWu/GCd4YYyJUq0gwHZNj+N9zBx6xr6i0jM05+SzZspfPN+3l82/38ObyHdwxZyXjMtpz3pAujO6dhsdT6c6mfD4y66psjDE1crtkshdoz5FLJkf08o7RPi99OiTSp0MiPz2xK6rKV9/v5/Uvt/Hm8u28uXwHfTskcsPpvZgwoIOTaMrnIzuU07TBG2NMBHCzZPINwJ3ALo5cMrlFrQcjIgw5vg1Djm/Dbydl8NaK7Tzy4Qau//uX9ExP4DcT+jK2TzsQj81HZowxLri5g/k50EdV94Q7mObC7/Nw7uAunDOoM/9ZuYOH5q3nyueWcEb/9vxfbCpeqyIzxpgauUkw3wO54Q6kOfJ6hEmZnRif0YGnF3zLXz5YxwZPLLHfb+E4VXc9z4wxppUKNVXMTcGnm4D5IvIWR64H80CYY2s2/D4P1/7gBCZldiT3yXbk79rK/774JX++IJOk6nqcGWNMKxdqJH9icPsOeB/wV9qXEP7Qmp/j2sbRv1dPesblM2/1Ls7+6wJWbW+VN3fGGFOjahOMqt6lqncB35Q/r7RvdeOF2LxIQjqJpfuY/bMTKSoJMOX/FvLeqp1NHZYxxjQ7bqbr/43Lfa1DfBqUFjCsUzRv3ngKfTsmce1LX/LyktY5/6cxxlQnVBvMBOAsoLOIPFzpUBJQGu7Amq3ywZaHdtOubQ/+ftWJTH9xKbe8spzcghKuOrVH08ZnjDHNRKg7mO3AEqAQWFppmwucEf7QmqnywZZ5zrow8dE+/nbZMCYO7Mi9b63m0Y82NGFwxhjTfISaTflr4GsReUlVW+8dy9Hi2zmPlQZbRvu8PHzxYKK8wn3vriU9MZoLhh3XRAEaY0zzEKqK7F+q+mPgKxE5ZgEUVW1RI/ldq2Y+Mq9H+PP5g8jJK+Y3r60gPSmG0b3TmiBAY4xpHkJVkf08+DgJOLuKrXWqmFH52KWT/T4Pj10yhN7tE7n2xaWs2GpdmI0xrVeobso7gk/HAn5V3VJ5a5zwmiFvlLMWTDXTxSTGRPHstOG0ifNz1fNfsCevqMrzjDGmpXPTTbkb8ISIbBSRf4nIDSKSFea4mrf49CrvYMqlJ8Xw1NRh7Msv4Zf/+ppAwJZYNsa0PjUmGFX9naqeDgwAFgC34PQma70SQicYgIxOSfxuUgafrMvm8U82NlJgxhjTfNSYYETkDhH5D/Ae0BO4Geji5uIicqaIrBWRDSJyWxXHRUQeDh5fLiJDaiorIvcEz10mIu+JSCc3sTSo+DRXi4799MTjmZjZkf/33jq+2Ly3EQIzxpjmw00V2RQgFZgHvAbMrdQ+U63gImWPAhOADOBiEck46rQJQK/gdjXwmIuy96lqpqpmAW8Cv3PxGRpWfFqNdzDgrDHzxykD6dImlhv+/hW5+SWNEJwxxjQPbqrIhuA09C8GxgErRGSBi2uPADao6iZVLQZmA5OPOmcy8Lw6FgEpItIxVFlVPVCpfDzO4meNKyENig5ASWGNpybGRPHIxUPIzivij++saYTgjDGmeXBTRTYAuAS4DLgQ2Ap86OLanYAMWAAAACAASURBVHHWkim3NbjPzTkhy4rI70Xke+CnVHMHIyJXi8gSEVmSnV3z3UatxB+eLsaNgV2SmXZSN/6x+DuWbrGqMmNM6+CmiuxPOFP0Pwz0U9UxquqmWqqq1biOvtuo7pyQZVX1dlU9DngJmFHVm6vqk6o6TFWHpaU18IDHhCOni3Hjl+N60zE5httfX0lJWaDmAsYYE+HcVJFNVNU/q+pCVa1NI8JWoPJ8KV1w5jdzc46bsgB/B86rRUwNo+IOxn2CiY/2MfOc/qzZeZBZC74NU2DGGNN8uLmDqasvgF4i0l1E/MBFOBNlVjYXmBrsTTYSyA12IKi2rIj0qlT+HKDxGzYSykfzu6siK3dG/w78sF97Hpq3nq378sMQmDHGNB9hSzDBCTJnAO/iLFD2L1VdJSLTRWR68LS3cZZk3gA8BVwXqmywzB9FZKWILAfGc3hKm8ZTPl2Mi67KR7trcn8A/vzO2oaMyBhjmp1qJ7ssJyLdVHXzUfuGq+oXNZVV1bdxkkjlfY9Xeq7A9W7LBvc3fpXY0aJiwZ9Yqyqycp1TYpl6Ulee/GQTN47tRc/0Vrn6tDGmFXBzB/OaiFTuwTUamBW+kCJEgrvBllW5+tQexPi8/PXD9Q0clDHGNB9uEsw1wBwR6SAiZwF/wVnpsnWrYT6yUFITopl6Ulfmfr2dDbvzGjgwY4xpHtz0IvsCuBFnqpiZwDhVtQXo63EHA4fvYh6xuxhjTAtVbYIRkTdEZK6IzAV+A8QBRcDTwX2tWz3uYMDuYowxLV+oRv77Gy2KSJSQDgV7obQYfP46XeLqU3vw/MItPPLheh66aHADB2iMMU0r1IJjH6vqx8AS4NPg8x1AMrCwkeJrvtqe4DzmrKvzJVITorl4xPG8uXwHew8VN1BgxhjTPLhp5P8EiAn2JPsAmAY8G86gIkKHgc7jrpX1usz5Q7tQGlDeXF7VRAXGGBO53CQYUdV8nGn7/6qq5wL9wxtWBEjtCd5o2LmiXpfJ6JRE3w6JvPbltgYKzBhjmgdXCUZERuHMXPxWcJ83fCFFCK8P0vvV+w4G4NzBnVn2/X42ZVtjvzGm5XCTYH6B04vs9eBULz2Aj8IbVoToMAB2rgSt35I0k7M6IwJzvrK7GGNMy+FmHMzHqnqOqv4p+HqTqt4Y/tAiQIdMyM+Bgzvrd5nkGE4+oR2vL9uG1jNZGWNMcxFqHMxDwceK8TCVt8YLsRlrP8B5bKBqsu/3FrBky756X8sYY5qDUONgXgg+2niY6rQP9nXYuQJ6javXpc4c0IE75qzktS+3Mbxb2wYIzhhjmla1CUZVlwYfP268cCJMbAokH98gdzDx0T7O6N+et5Zv586zM4iJsn4UxpjIVm2CEZEVHLvEcQVVzQxLRJGmvKG/AZw1sCNzlm1nxbZcu4sxxkS8UFVkk4KP5eu1lFeZ/RSw5RjLdRgI696BkgJnnZh66NzGKZ9zsKghIjPGmCYVaqqYLaq6BThZVX+tqiuC223AGY0XYjPXfgBoAHZ/U+9LtUuIBmCPTRtjjGkB3IyDiReRU8pfiMhJQHz4QoowHYI9yRqgmqxNnDNp5p48SzDGmMhX45LJwJXALBFJxmmTyQWuCGtUkSSlG/gTGqSh3+/zkBTjY+8hqyIzxkS+GhNMsDfZIBFJwpmXLDf8YUUQj8epJqvnnGTlUhOiybEqMmNMC+CmigwAVT1gyaUaHQbArlX1njIGIDXez16rIjPGtACuE4wJof0AKDoA+7fU+1KpCX72WBWZMaYFsATTEMrXhmmAhv628dG2+JgxpkVw08iPiAwAMoCY8n2q+ny4goo46f0AcRr6+02q8fRQ2iX42XuomLKA4vVIw8RnjDFNoMYEIyJ3Aj/ASTBvAxOABYAlmHL+eGcBsgZo6G8b7yegsD+/mNTguBhjjIlEbqrIzgfGAjtVdRowCLDffEfr0DA9ycqTilWTGWMinZsEU6CqAaA02FV5N9AjvGFFoPYDnEb+wgP1ukxqvDPYMsd6khljIpybBLNERFKAp4ClwJfA4rBGFYnKG/p3rarXZVITnARjdzDGmEjnZqDldcGnj4vIO0CSqi4Pb1gRqPLiY11H1fkyqfHl85FZV2VjTGRz1YusnKpuDlMckS+pE8S2hZ31y71t4qIAm4/MGBP5bBxMQxFpkLVhfF4PbeKi7A7GGBPxLME0pPYDYfdqCJTV6zJt4/3WBmOMiXhuB1q2AY6rfL6qfhmuoCJWhwFQWgB7NkJa7zpfJjUh2nqRGWMinpuBlvcAlwMbObyEsgKnhy+sCFXRk2xF/RJMvJ/1u/MaKChjjGkabu5gfgycoKr2J3VN2vUBT5Qz4HLAeXW+TGqCn0WbrA3GGBPZ3LTBrARS6nJxETlTRNaKyAYRua2K4yIiDwePLxeRITWVFZH7RGRN8PzXg2N0mgefH9L61Luhv218NPvySygtCzRQYMYY0/jcJJg/AF+JyLsiMrd8q6mQiHiBR3HmLssALhaRjKNOmwD0Cm5XA4+5KPs+MEBVM4F1wG9cfIbG035AvVe3bBccbLkvv6QhIjLGmCbhporsOeBPwAqgNn9SjwA2qOomABGZDUwGvql0zmTgeVVVYJGIpIhIR6BbdWVV9b1K5RfhzJXWfHQYAMtnw6EciG9Xp0tUHmyZlmjTvhljIpObBJOjqg/X4dqdge8rvd4KnOjinM4uywJcAfyzqjcXkatx7oo4/vjjaxN3/VSsDbMCThhTp0u0Dc5HZitbGmMimZsqsqUi8gcRGSUiQ8o3F+WqWszk6DWFqzunxrIicjtQCrxU1Zur6pOqOkxVh6WlpbkIt4G0L+9JVvdqsvIqshwbC2OMiWBu7mAGBx9HVtrnppvyVpyxM+W6ANtdnuMPVVZELgMmAWOD1WvNR3wqJHasV0P/4TsY60lmjIlcbia7rFs9D3wB9BKR7sA24CLgJ0edMxeYEWxjORHIVdUdIpJdXVkRORO4FRitqvl1jC286tnQnxLnxyOwx+5gjDERzM1Ay2TgTuC04K6PgbtVNTdUOVUtFZEZwLuAF5ilqqtEZHrw+OM4K2SeBWwA8oFpocoGL/0IzoJn74sIwCJVne7+IzeCDgNg03woLXa6LteS1yO0ifNbgjHGRDQ3VWSzcMbC/Dj4+lLgGWBKTQVV9W2cJFJ53+OVnitwvduywf09XcTctDoMhEAJZK+Bjpl1ukRqgp89VkVmjIlgbhLMCapaeVj6XSKyLFwBtQiVG/rrmGDaxvttyn5jTERztWSyiJxS/kJETgYKwhdSC5B6AsQkw7ef1v0SCdE2o7IxJqK5uYOZDjwfbIsB2AdcFr6QWgCPF/qcBWvfhrIS8EbV+hLt4v3kWBWZMSaChbyDCU7ZcomqDgIygUxVHWxLJrvQ7xwo3A/fflKn4m3jozlQWEpxqc1HZoyJTCETjKqWAUODzw+o6oFGiaolOGEMRMXD6hqnbatSasV8ZFZNZoyJTG7aYL4KTnB5qYhMKd/CHlmki4qF3uNhzVt1WuEyNTjY0hr6jTGRyk2CaQvswRm5f3ZwmxTOoFqMfufAoWz4blGti6YmHJ7w0hhjIlG1jfwi8idVvRV4W1VfbsSYWo5e48Eb7VSTdTu5VkUrpouxnmTGmAgV6g7mLBGJormttxJJohOg51hY/QYEatdYXzHhpVWRGWMiVKgE8w6QA2SKyIFK20ERscZ+t/qdAwe2wfavalUsKSYKn0fYa1VkxpgIVW2CUdVbVDUZeEtVkyptiaqa1IgxRrY+Z4LHB6v/XatiHo/QITmG7/bamFZjTGSqsZFfVSc3RiAtVmwb6D4avpkLtVxZoG+HJNbssJtFY0xkctOLzNRXv7Nh37e1nsK/X8dENuUcorCk9t2cjTGmqVmCaQx9J4F4nMb+WujXMYmygLJhd16YAjPGmPCpMcGIyM/d7DMhJKTB8Sc51WS10LdDIgCrrZrMGBOB3NzBVDWx5eUNHEfLl3EOZK+GnPWui3RNjScmysPqHQfDGJgxxoRHtQlGRC4WkTeA7sGpYsq3j3BG9pva6Buc/OAb973JvB6hT4ck1uy0OxhjTOQJNV3/QmAH0A74f5X2HwRsNuXaSu4MnYc5o/pPu9l1sX4dEnl31U5UleAS0cYYExFCjYPZoqrzVXWUqn5caftSVUsbM8gWI+Mc2PE17Nviuki/jknsyy9h90EbcGmMiSxuGvmniMh6Ecm1kfz11O9s57EWvcmsod8YE6ncNPL/GThHVZNtJH89te0B7QfWao2Yvh2cr9oa+o0xkcZNgtmlqqvDHklrkXEOfP85HNjh6vTkuCg6p8RaQ78xJuK4STBLROSfwV5ltuBYffU7x3lc86brIn07JLLG7mCMMRHGTYJJAvKB8diCY/WX3hfa9a5Vd+W+HRPZmJ1HUalNGWOMiRyhuikDoKrTGiOQViXzQvjwHlg1B/r/qMbT+3VMojQ4ZUz/TsmNEKAxxtSfm15kvUXkAxFZGXydKSJ3hD+0FuykG6HzUJh7I+zbXOPp1tBvjIlEbqrInsJZ1bIEQFWXAxeFM6gWz+eH854GFF65EspKQp7eLTWOaJ/Hpu43xkQUNwkmTlUXH7XPBlrWV9vucM7DsG2JU10Wgs/roU+HRNbstDsYY0zkcJNgckTkBEABROR8nClkTH31PxeGXg7//QtsmBfy1H4dkli+db+tDWOMiRhuEsz1wBNAXxHZBvwCuDasUbUmZ/4R0jPgtWvg4M5qTzt7UCcOFJbyxtfbGzE4Y4ypOzdLJm9S1R8CaUBfVT1FVTeHPbLWIioWzn8Gig/Baz+DQNV3KCf3TKV3+wSe+e9mtJZLLxtjTFNw04ssRURuBO4Bfi8iD4vIw+EPrRVJ7wtn/Rm+/QQWPFjlKSLC5Sd155sdB/hi875GDtAYY2rPTRXZ20A3YAWwtNJmGtLgS2HAefDR/8J3i6o85dzBnUmJi+KZ/37byMEZY0zt1TjQEohR1ZvCHklrJwKTHoJtS52uy9M/hbi2R5wS6/dy0fDjefKTjWzdl0+XNnFNFKwxxtTMzR3MCyLyMxHpKCJty7ewR9YaxSQ57TF5u2DuDVBFW8vUUV0REV74zP2aMsYY0xTcJJhi4D7gMw5Xjy0JZ1CtWuch8MOZzmSYi5865nCnlFjO7N+Bfyz+jvxiG45kjGm+3CSYm4CeqtpNVbsHtx5uLi4iZ4rIWhHZICK3VXFcgp0GNojIchEZUlNZEblARFaJSEBEhrmJI+KMuh56nQHv3Q47jl2d+opTunOgsJSnP7W2GGNM8+UmwazCmU25VkTECzwKTAAygItFJOOo0yYAvYLb1cBjLsquBKYAn9Q2poghAj96DOJS4ZVpUJR3xOGhXdswYUAHHp2/gW37C5ooSGOMCc1NgikDlonIE+VdlF12Ux4BbAiOoykGZgOTjzpnMvC8OhYBKSLSMVRZVV2tqmtdfr7IFZ8K5/0N9m6Ct28+5vDtE/sB8L9v2VpwxpjmyU2CmQP8HlhI7bopdwa+r/R6a3Cfm3PclA1JRK4WkSUisiQ7O7s2RZuPbqfAab+Gr/8By/5xxKEubeK4dnRP3lqxg4UbcpooQGOMqZ6bkfzPAf8CFqnqc+Wbi2tLVZdzeY6bsiGp6pOqOkxVh6WlpdWmaPMy+tfQ9RR461eQs/6IQ9eM7kGXNrHMfGMVJWWBJgrQGGOq5mYk/9nAMuCd4OssEZnr4tpbgeMqve4CHD2RVnXnuCnbOni8cN5T4IuGl6dBSWHFoZgoL3dMzGDdrjyeW7i56WI0xpgquKkim4nTJrIfQFWXAd1dlPsC6CUi3UXEj7OGzNGJaS4wNdibbCSQq6o7XJZtPZI6wbmPw64V8P5vjzh0Rv/2/LBfOn9+Zy0rt+U2UYDGGHMsNwmmVFWP/s1VY3WVqpYCM4B3gdXAv1R1lYhMF5HpwdPeBjYBG3AWNrsuVFkAETlXRLYCo4C3RORdF58h8vU+A0ZeD4ufhNVvVuwWEf58/iDaxvuZ8fcvOVgYevEyY4xpLFLTzLwi8jTwAXAbcB5wIxClqtNDFmxGhg0bpkuWtICxoaVF8PR4Z5nl6Qsg5XAt4heb93LRk4s4a2BHHr4oC5GqmrGMMcY9EVmqqnUeb+jmDuYGoD9QBPwdyMVZE8Y0Nl80nD/LmdL/1Suh7PBI/uHd2nLTuN688fV2/rH4+xAXMcaYxuGmF1m+qt6uqsOD2x2qWlhTORMmqSfA2Q/B95/D/P894tC1o0/gtN5pzJy7ioUbreuyMaZpuelF9r6IpFR63abVtHs0VwPPd6b3//QB2PhRxW6PR/jrRYPp1i6Oq59fao3+xpgm5aaKrJ2q7i9/oar7gPTwhWRcmfAnaNcbXrsa8nZX7E6Oi+K5K0aQFOPj8me+YMueQ00YpDGmNXOTYAIicnz5CxHpSi0HPZow8MfDBc9A0QF4/RoIHB5o2TE5luevPJGyQIBLn17MrgNWo2mMaXxuEsztwAIReUFEXsCZZPI34Q3LuNK+P5z5B9j4ISw8cnq4nukJzLp8OHvyijjvsYV8m2N3MsaYxuWmkf8dYAjwT5wpY4aqqrXBNBdDp0HGZPjwHvj+iyMODT6+Df+4eiT5xWWc/9hCa5MxxjQqN3cw4MyovBuni3KGiJwWvpBMrYjA2Q87o/1fuQIK9h9xOLNLCq9MH0VMlJeLnlzEgvXWu8wY0zjc9CK7Cqda7F3gruDjzPCGZWolNsVZavng9iqXWu6RlsCr155Ep5QYps76nEc/2kAgYM1oxpjwcnMH83NgOLBFVccAg4EInf++BesyDMb+DlbPhVlnwJaFRxzukBzDa9edzMTMTtz37lquen4J+/OLmyhYY0xr4CbBFJYPrBSRaFVdA/QJb1imTk66Ec7+C+zbAs9MgJd+DLtWVRxOiPbx8EVZ3D25P5+uz+asv3zKJ+vsbwVjTHi4STBbgwMt5wDvi8i/aa1T5zd3IjD0crjxKxh7J3y3CB47GV67xkk6OJNjTh3VjVemn0SM38vUWYu55eWvyc23STKNMQ2rxskujzhZZDSQDLwTXMo4IrSYyS5rK38vLHgQPn8CUBh+FZx6s7McM1BYUsbDH6zniU82kRrv5zdn9WXyoM54PDZRpjGm/pNdVptgRKRtqIKqureub9rYWm2CKZe7Feb/AZb9HaLi4eQbYeR1EJ0AwMptudz22nJWbjtAZpdk7piYwYjuIf/5jTGtQDgTzLeEWL5YVXvU9U0bW6tPMOV2r3HGy6x5E+LTneWYh1wGPj+BgPL6V9u479217DxQyNi+6dw4theDjkup+brGmBYpbAmmJbEEc5TvF8O8mbDlv9CmG5z+W+g/BTweCorLeHrBJp769FtyC0o4rXcaM8b0ZHi3NrbGjDGtTKMkGBGZApyCc0fzqarOqesbNgVLMFVQhfXvO4lm9yrokOkkmhNOB6+PvKJSXvhsC3/7dBN7DhUzoHMSl43qxtmDOhET5W3q6I0xjSDsCUZE/g/oCfwjuOtCYKOqXl/XN21slmBCCJTBipfhw99D7ncQkwInjIGe46DnWPKj2/Hal9t4buFm1u/Oo228n/OGdOaCYcfRu31iU0dvjAmjxkgwq4ABGjxRRDzAClXtX9c3bWyWYFwoLYK1b8P6ebBhHuTtdPZ3GAg9x6E9x7Ko+ASe+3w781bvojSgDOqSzHlDuzBhQEfSEqObNn5jTINrjATzGvBLVd0SfN0V+KOqXlzXN21slmBqSRV2rnASzYZ5zngaLYPoJOjxA/KO+wFv5PfnuZXFrNl5EI/AyB6pTMrsxA/7pZOeFNPUn8AY0wAaI8F8jDNVzOLgruHAZ0A+gKqeU9c3byyWYOqpMBc2fQwb3nfucA4Gx9mm92dvp1P5oCSTpzans26PMzRqUJdkxvZrz+jeaQzonIzXxtUYE5EaI8GMDnVcVT+u65s3FkswDUgVdq8OJpv3nbubQAnqTyCv00ksjRrGi3t6M2+7H4CUuChOOiGVUT1SGd69Lb3TE20gpzERojESTDxQoKoBEekN9AX+o6oRM7eIJZgwKjoI337iJJsN8yD3ewDKko4jJ7Yba0o78t/cVJYcas8G7QQxKQzp2oZBXVIYdFwyAzunWPuNMc1UYySYpcCpQBtgEbAEyFfVn9b1TRubJZhGogo565xks/1LyF4He9ZD6eElmw/42rJJO7OiuAPrAp3ZoJ3ZH9ed9I7H069TMn07JNKrfQInpCVYd2hjmlh9E4zPzXuoar6IXAn8VVX/LCLL6vqGpgUTgbQ+zlYuUAb7tzjJJmctSdlrycpey6CcRUjRQeecUji4NZ71WzqxPtCZ17Uzm+hEflJP4tO70z0tkR5pCXRrF0fX1Hg6JsVYNZsxEcBVghGRUcBPgSuD++xPS+OOxwttezhbnzMrdosqHNwJ2WsgZx2J2WvJ2r2GzOyV+ArmOycVQOEWP7s2tyFbk8nRZD7SJPZJCqVxaXgS2xOT0oGE1E4kp3WiY1o7OqXEkp4YYx0LjGkG3CSYXwC/AV5X1VUi0gP4KLxhmRZPBJI6OtsJYwBn7QgPOLNAZ6+FnLXE5Kzn+IM7ab9/J2UHd+HNX09MyX4owtlygA3OJfM1mhxN4mtSOOhrQ6E/ldK4NIhPw5fUgZiU9sS17URyWhfatW1LSpzf7oSMCSObi8xEnrISOJQDh3ZDXjYF+3eQt2c7Bft2UnZgJ578bPyFOcSV7CUxcAAPx/6MF6ifHJLZ70khz9eWQn9bimPaEYhPQxLaE5XUHn9KB+LbdiQpuS1t4qNJifPbnZFpVcLWBiMiD6nqL0TkDTj2f2gkjH8xLZQ36vDdDxAb3KpUVgr5eyjYt4Pc7G0c2redov07KTuwG0/+bqIL9pBSnE1C/jqSDuXi3RM45hJFGkU2yXyvSeR6Ujjoa0txVDLqT0CjE/FEJ+CNTSIqNgl/XDLR8cnEJaQQl5hCXGIySfFxJMb48HndrO9nTMsRqorsheDj/Y0RiDFh4fVBYntiE9sTe3xW6HMDZZC/h6L9O8nbs538fTsozt1J4OBu5NBuUgpyaF+UQ1zxZmIL8/AXuFtzr1Cj2E8s+cRS4ImjyBNLkTeeUl88Zb54yqISUH8iRCciMQl4Y5LwxiYekbBiEpKJTUgmPi6eOL/XZrY2EaHaBKOqS4OPH4tIWvC5LeBuWi6PFxLSiU5IJ7pLJqk1nV9W4owDKs4jUHCAgkO55OftpzAvl6JDuZQWHKCs4ACBwoNQdBApzsNTcoiY0jwSS/fhL9pGbEE+sVpALEWuQixSH/uIIZ9Y8iWOQo+zlXjjKfXFEvBGo74Y8MVCVAzii0X8sXiiYvFEx+Lzx+H1x+KNjiMqOg5/TBz+mHj8sfH4Y+KIiY0nJiYGn0csiZl6C1VFJsCdwAycRcc8IlKK01X57kaKz5jmyxsFcW0hri2eFIjH2eqkrJSyooMU5OVScHAfhYdyKco/SGl+LiUFBygrOIgWHYCiPKToIN6SQ3hKy5NVHtFlu/CXFODXYvxaRAx1X9G8VD0cwk8RforFT7FEUyx+SiSaEk80pZ5oyjzRlHljKPPGVCS1gDfGSW5RseCLQaKczeuLxhMVjdcXhTcqGm9UNL6oaKKi/Pj8wed+53mUP4Yovx+/Pwa/z2edMCJcqCqyXwAnA8NV9VuAYA+yx0Tkl6r6YGMEaEyr4PXhjWtDQlwbEtK71f96qlBahJYUUFx0iOL8fAoL8yguLKCk8BAlRYcoLSqgrKiAsuJ8ykoKCRQXQElwKy1EgpunrBBvWRHeQCGxgSJ8ZYeIKi0iSouCCa2YaIqJorT+cVdSqh6K8FGKlxLxUUoUZXgplSjKxEepRBEQL2USRcATRUB8zqMnCq20BTxRzh8D3ijU60cqXvvB58cTfO7x+fH4ohCvH/H58UZFB/f58fr8eKKcRydJ+vFFReP1+fFFRxPli8bn9+PzWvVlZaESzFRgnKrmlO9Q1U0icgnwHmAJxpjmSsSpIouKITquDdFtIOyr9wTKoKSAQHEBxUX5wURWQElxEaUlRZSVFB/xGCgtprSkmEBpMVoafCwrhrJitLQEDZRAqfNaAiUQKMUTfC6BUjyBEjxagif43K/5eLUUr5bi01J8WoKXMnyUEKWlRFGKjzJ8cmxHjoZSol5KypMiPkrFR2mlxzLxVSTHMvEREB9lEkWZRKEeL4gXDW54Dj8SPIbH5zx6fSAe8PiQ4HHx+hDxIV6f89rjnO8pP+ZxjpW/rnj0+hCvF4/Xhyd4jtfrw+tzM4oltGq7KYvISlUdUNtjzZGIHATWNnUczUQ7nNEjxr6Lyuy7OMy+i8P6qGqd/zYJlaJCVeLWvYK3aaytT1/ulkRElth34bDv4jD7Lg6z7+IwEanXAMJQCWaQiByo6j0BW1HKGGNMSKG6Kdt8Y8YYY+qstQwtfrKpA2hG7Ls4zL6Lw+y7OMy+i8Pq9V20irnIjDHGNL7WcgdjjDGmkVmCMcYYExYtOsGIyJkislZENojIbU0dT2MSkeNE5CMRWS0iq0Tk58H9bUXkfRFZH3xs09SxNhYR8YrIVyLyZvB1q/wuRCRFRF4RkTXBn49Rrfi7+GXw/8dKEfmHiMS0lu9CRGaJyG4RWVlpX7WfXUR+E/xdulZEznDzHi02wYiIF3gUmABkABeLSEbTRtWoSoFfqWo/YCRwffDz3wZ8oKq9gA+Cr1uLnwOrK71urd/FX4B3VLUvMAjnO2l134WIdAZuBIYFB457gYtoPd/Fs8CZR+2r8rMHf3dcBPQPlvm/4O/YkFpsggFGABtUdZOqFgOzgclNHFOjUdUdqvpl8PlBnF8inXG+g+eCpz0H/KhpImxcItIFmAj8rdLuVvddiEgScBrwNICqFqvqflrhdxHkA2JFxAfEAdtpJd+Fqn4C7D1qd3WffTIwW1WLgnNTbsD5laWK1AAABwRJREFUHRtSS04wnYHvK73eGtzX6ohIN2Aw8DnQXlV3gJOEgPSmi6xRPQT8Gqg8EVVr/C56ANnAM8Hqwr+JSDyt8LtQ1W046119B+wAclX1PVrhd1FJdZ+9Tr9PW3KCqWpK01bXJ1tEEoBXgV+oalUzM7R4IjIJ2F2+xlEr5+P/t3euMXZNURz//du0oxj9gEgbj1bVI4IhGhQxpXxohEYqRElDPEIolfrg8aESkpakBFHxaBo0BJ3UeBRJvZ990BoUiVZiPNpGvFPKdPmw1nB6e2faO3rn6r3rl5xk3/04e+1177nrrL3PWRuOAuaY2ZHAb9TvFFCvxPrCmcBIYDiwSwTzTbakT/+n9WxgOoF9Cp/3xt3fhkHSINy4zDeztsheK2lYlA8D1tVKvn7keOAMSV/iU6UnS3qUxtRFJ9BpZu/F56dwg9OIuhgPrDGz9Wb2J9AGjKUxddFNT2Pv0/9pPRuYpcBoSSMlDcYXqNprLFO/ERvGPQSsMrPZhaJ2YEqkpwBP97ds/Y2ZXW9me5vZCPx38LKZnU9j6uI74CtJB0XWKcAnNKAu8KmxYyXtHNfLKfhaZSPqopuext4OnCupSdJIYDSwZGsnq+s3+SVNwOfeBwJzzezWGovUb0g6AXgD6ODfdYcb8HWYJ4B98QvsbDMrXeirWyS1AtPN7HRJu9OAupDUgj/sMBhYDVyI32w2oi5uBs7Bn7r8ALgY2JUG0IWkx4BWfHuCtfgOxgvpYeySbgQuwnV1jZkt2mof9WxgkiRJktpRz1NkSZIkSQ1JA5MkSZJUhTQwSZIkSVVIA5MkSZJUhTQwSZIkSVVIA5PUBEldklZEJNuVkq6VNCDKjpZ0Vy9tR0g6r/+k3aL/qRGFeH6tZOgrkiZuj6Cvko6U9GCkZ0ia3sfzDJb0esQCS+qMNDBJrdhgZi1mdihwKjABfw4fM1tmZlN7aTsCqJmBAa4AJpjZ5GLmDvInORGPLr7N9DCuG4C7/6swEYh2Mf4uSlJnpIFJao6ZrQMuBa6U01rYs+Wk8HRWRHDGZmAmcGLkTQuP5g1J78cxNtq2Snq1sPfJ/HhjG0ljJL0d3tMSSc3y/WJul7RU0oeSLiuVVdJ9eMDI9uh7hqT7Jb0EPCxpP0mLo/1iSftGu3mS5sj36Fkd45obntC8cnqpRMYY62uSnpD0uaSZkiZHuw5Jo0IvZwC3h+5GxfGCpOWhw4ML8s6W9Aowq0SuZuBwM1tZRuZLJC2SNCR0f0d4KKtiPG3yvUZuKTRbCEwuPVdSB5hZHnn0+wH8WibvB2Av/O3iZyPvGeD4SO+KB2v8pzzydwZ2ivRoYFmkW4Gf8LhJA4B3gBP49w32MVFvtzjvpcBNkdcELANGlpHzS2CPSM8AlgNDCvJOifRFwMJIz8PjoAkPsPgzcFjItRxoKemjIhljrD8CwyL/a+DmqHc1cGdBjkmFfhYDoyN9DB5Gp7ves8DAMuMfBywofJ4BTAeuxEOKNEX+q8CsggzfFOTrBHaPsoHA+lr/JvPY/seO4NInjUO5iK1vAbNjvaPNzDrDCSkyCLgnQqB0AQcWypaYWSeApBX49NpPwLdmthTAIsq0pNOAwyVNirZDcYO1Zityt5vZhkgfB5wV6UeA2wr1njEzk9QBrDWzjuj345BrRaHuQRXKuBFYahFqXdIXwEtRpwM3Cpshj7Q9FniyoNOmQpUnzayrzHiH4SH/i1yAG42J5oEju+mO/9cBfFyQbzUePPF7M+uStFFSs/neRUmdkAYm+V8gaX/cOKwDDunON7OZkp7D12jelTS+TPNpeCylI3CP4PdC2R+FdBf+mxflQ40LuMrMXqxQ/N96KSv20y3LphK5NrHltViRjPIYa6XnLPZX7lofAPxoZi09yN7TuDYAO5XkfQS04N5i0SBv65ib2Px7S+qAXINJao6kPYH7gHvMzErKRplZh5nNwqeDDgZ+AZoL1Ybid/ub8DvprW3l+ikwXNKY6KM5FrJfBC6Xb3OApAPlm3FVwtt4xGbwdYU3K2xfbRn/0V14RWsknR3nkqQjtuEcq4ADSvI+AC7D16aGVyAP8qCj60s8n6QOSA8mqRVDYspqEB6d9RFgdpl610gah3sfnwCL8LvfvyStxNcK7gUWxB/lK/TuUWBmGyWdA9wtaQh+Rz4ejzA8Ang/HgZYT+Xb5U4F5kq6LtpfWGH7asv4OPCApKnAJNwIzpF0E/5dPA5ssXhfItunkoaWTmmZ2Zvyx5Wfk3RqBTKNA56voH6yg5DRlJMkqRhJ04BfzOzB7XCuNuB6M/vsv0uW/J/IKbIkSfrCHDZfU+kT8s0AF6ZxqU/Sg0mSJEmqQnowSZIkSVVIA5MkSZJUhTQwSZIkSVVIA5MkSZJUhTQwSZIkSVX4GziVNXgUY0s9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(dg_cm,np.abs(ug_cm),label='Disk load response')\n",
    "plt.plot(km_ld, np.abs(u_m_kg*(area_disk*1e6)*rho_w),label='Theoretical response, only valid in the far field')\n",
    "plt.xlim([0,100])\n",
    "plt.ylim([1e-4,0.005])\n",
    "plt.legend()\n",
    "plt.ylabel('Displacement from a disk with 1m water (m)');\n",
    "plt.xlabel('Distance from cemter (km)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#If you only want the unit load Greens's functions (e.g., for ice sheet applications), stop here and save the data using \n",
    "# sio.savemat('data/greens_functions_025_deg.mat',{'dg_cm':dg_cm, 'ug_cm':ug_cm})\n",
    "\n",
    "#Otherwise the next step is to set up a Greens function matrix (G) so that G.m = d,\n",
    "#where m is the load and d is the displacement. For terrestrial hydrology applications, \n",
    "#I need a good land/ocean mask – I use GMT here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#First, define the domain of interest\n",
    "min_lat = 52\n",
    "max_lat = 72\n",
    "min_lon = -175\n",
    "max_lon = -132\n",
    "\n",
    "gmt_R = '-R' + str(min_lon) + '/' + str(max_lon) + '/' + str(min_lat) + '/' + str(max_lat)\n",
    "os.system('gmt grdlandmask -I' + str(grid_cell) + ' -r '+ gmt_R + ' -N0/1/1/1/1 -Gdomain_mask.nc');\n",
    "os.system('gmt grdmath ' + gmt_R + ' AREA domain_mask.nc MUL = domain_mask_area.nc');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from netCDF4 import Dataset, num2date\n",
    "domain_mask_area = Dataset('domain_mask_area.nc')\n",
    "\n",
    "lat_grid = domain_mask_area.variables['lat'][:]\n",
    "lon_grid = domain_mask_area.variables['lon'][:]\n",
    "area_grid = domain_mask_area.variables['z'][:]\n",
    "lon_grid,lat_grid = np.meshgrid(lon_grid,lat_grid);\n",
    "\n",
    "glon = lon_grid[area_grid!=0]; glat = lat_grid[area_grid!=0];\n",
    "area = area_grid[area_grid!=0];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydeZwU1bX4v6eqe1YYdlR2EMQFFRVRgxoXjEiIGjUu0bglMSbG+J5JDL7kJVF/vqcv0aiJWXCLURH3iERxRXEJCIgiq8i+7+ssvVSd3x9VPdPT091T3dPT0z1T3w+XqeXWrVPVVXXuPffcc0VV8fHx8fHpmBhtLYCPj4+PT9vhKwEfHx+fDoyvBHx8fHw6ML4S8PHx8enA+ErAx8fHpwPjKwEfHx+fDkygrQXw8fHxaU+IiAnMBTao6gQReQYY7u7uCuxW1ZEiMghYAixz981S1evzLa+vBHx8fHxyy004H/cqAFW9JLZDRO4B9sTlXaGqI/MrXmN8c5CPj49PjhCRfsDXgYeT7BPgYuDpfMuVjoJsCfTs2VMHDRrU1mL4+PgUOPPmzduuqr1aUsbZp1fqjp2Wt/MtCL2uquPSZLkPuAXonGTfKcAWVV0et22wiMwH9gK/UtX3PYqdMwpSCQwaNIi5c+e2tRg+Pj4FjoisaWkZO3ZafPz6AE95zYOWHyoi8R+nSao6yZVlArBVVeeJyGlJDr+Mxq2ATcAAVd0hIscB/xSRI1R1b1YXkiUFqQR8fHx88oUCNrbX7NtVdVSKfWOAc0VkPFAGVInIk6p6hYgEgAuA4+rPqxoCQu7yPBFZARyC06mcN/w+AR8fnw6NokTU8pTSlqN6q6r2U9VBwKXAO6p6hbt7LLBUVdfH8otIL9eTCBEZAgwDVrbGNabDbwn4+Ph0eDJoCWTLpTTtED4VuF1EooAFXK+qO1tbkER8JeDj49OhURQrxyH1VfVd4N249auT5HkBeCGnJ84CXwn45JRwXZgv5q1k7459mAETMYRoOIphGJglJtFQFFWlsqqcwUcOpKpHMicKH5/8YtNx51XxlYBPznhnygf84ft/JVwbxrabeakEAsEA5/3obH5wz1U4LtQ+PvlHAasDKwG/Y9gnJ6xauJZ7v/sX6qpDzSsAAIVoOMq0SW/x6kNvtb6APj5psFFPqT3itwR8WoSqMveNBdz/o4cI1YYzPj5UE+KPNz7KtIfeYcSY4Zz/43H0HXogm1Zu4aU/TefzD5ZRs68GwWkpKOosC3TuWskJXz+WCdeNpWuvqlxfmk8HQYFIB55m11cCPi3ivh89wttPvU+oujbrMqxIlBWfrWHFZ2uY/ti7fOfXF/GP258n7EGpfPHJKl68/zX++OHt9B12UNYy+HRcFPXNQT4+8VhRi7qaEHXVIUJ1Yepq3eVad70mRF1NiMWzl/PO5A8I10VAWvAoxfUHhGrDPPLLpz0pAAC1leo91fz55n9kf36fjo2C5TG1R/yWgE89tfvrePCnT/D2lI+wreb9ptWyIJavJR27htm4XC99CvH5FebPWJT9+X06NM6I4Y6LrwR86vnNxfex8KNlnhRAIiKCmgGwbVCvxwuYZk48g0zTb9T6ZItg0XG903wl0MHZvGYb//7Xpyz+eDmff7QMO5qBAjCMhpYAjiLANAEz9TGthGUrf/zpk1jRKKZpYgRMouEoAhhBAztqo6qYQRO1HZOXYRoYAcPNJwSCJpab74ABPRg19kiGHuUtsJhP8eJ0DPtKwKcD8vTvpzH5d9OIhKOobTf6oHtBRNCACVFvYXg9k2mZIkRt5V+PvptTMZ76v1f4yteP4RcPfR/D8Fsa7RVnnEDHVQL+k91BWf7ZGqbc+yqRcNTZkKVJRgwDggHnw20aYIiT4pcNo2FbIICUBJ1jTKMhn2k4ZQQDKcpMUV4wgAQDrTLYLBKKMmv6Z7zz7Kycl+1TWNgqnlJ7xG8JtGM+fmsh/7h7Kuu+2IwVtUEUVcfP3gpH0bjatoigIk4va4aISMZKpMF0lKMy05TVEkI1YX5/w2Pc99OnQAB1ao6O0lHUjokpiECgJMAhIwdy9a3ncdjxQ1pFJp/c4rcEmkFEhovIp3Fpr4j8h4h0F5E3RWS5+7dbiuPHicgyEflSRCbm/hJ8kvH6Ux9y53cnsWLBOsJ1EayohRWxsaPq2L3tJKYf03Rq2cWGaTqth9ZCnXtmRWysqI0dtbEi7v20nG1W1CIasairDrHgwy+49Vv3MX/m0taTySdnKIKF4Sm1R5ptCajqMmAkgBv7egPwEjAReFtV73I/7hOBX8Qf6+Z/EDgLWA/MEZGpqro4p1fh04htG3fxt18/5/jvp8IwwG5sd6+vnbdSrbpoyWIMRKg2wv03P8kN/3cZVtSuD6ZnRaIgQqDEJBq2UFsJlpjYtqNojICBaRpEI1FUIVgaxIpa9DqwK4MO6+PHWGol2qupxwuZmoPOBFao6hoROQ84zd3+OE7Y1F8k5B8NfKmqKwFEZApwHuArgVbAilrcd8vTvPPCHOy6ZgZbiWtfz9Anv8Mhbt9DFmxZt4NfX/m3nIhhBgx69+3O7U9cT78hvXNSpo+DIoS141Z8Mn264ydGOEBVNwG4f5M9mX2BdXHr691tPq3A0w+8wcypn3hy8xQRJOB2vsZXgsRNpkEHNpM6H/+ACYHcjGNoKVbUZtOa7Uy8+I9YWYzj8EmNM1jM8JS8ICKmiMwXkWnu+m9FZEOcSX18XN5bXVP5MhE5u3WuMD2eWwIiUgKcC9yaQfnJ3p6kVU8RuQ64DmDAAN83OxteemgG4TrH3IDHTl4xjNT9AL5ZqGW0Qj/Fzi17ufKk31JWUVq/TUSo6l7JV752JOMuPYlOXcqbHLd53Q6mPfkhG1dt46gTh3LWt06gsnNZzuUD2LZpN/964kPWfrmZw44dxLjLTqJzl4r6/arKnBmLeW3yv9m4ZjtWxIncY5gGPQ7owhnnH8dp5x5LSVmwVeRLRo47hm8ClgDxUQ3/oKq/j88kIofjVKyPAPoAb4nIIarNzGOZYzIxB50DfKKqW9z1LSJykKpuEpGDgK1JjlkP9I9b7wdsTFa4qk4CJgGMGjXKt1FkyFMPvE7NvrqGDaYJlpWVt49PDkinXFuAqrJz854mZW9YtY1ln67lxYfe5YFpN9PzwK71+z6fvYJfXz2JqNt5/cn7y3h+0gz+NO2ndO2Z20l9vliwlomX/ZloJEokbDFv5jJeeOhd/vjKzfTq0w1V5Z6fTub9f31KOBRtcvz6FVtZNGclLz82k3tevImy8pKcypcMVcHS3PxWItIP+DpwJ3BzM9nPA6a4E86vEpEvcUzo/86JMB7J5Movo/EcmVOBq9zlq4CXkxwzBxgmIoPdlsSl7nE+OWTbpt0882BCTH4RCAScFOvsTZYMI/3+tkiFKFMmKXbPW8uMlKJc27LZs2s/j909jXAoQiQcIVQX5p6fTqauNkw04lQwQ7URdu/Yz1P3v54TccKhCJFIlFBdhD/8fAq11SEiYedc4boI+3ZV86gr05JPVvP+q58lVQAxohGL9au2MX1K/r6FNuIpeeA+4BaahiP6sYgsEJFH4zwpC8Jc7qklICIVOB4+P4jbfBfwrIh8F1gLfMvN2wd4WFXHq2pURH4MvI4TS+BRVfUjfeWYFx9+1xn0ZRpNR/1m4cPvU8AIaX9P21LeeWke70xb4GxQhWjTD64Vsfjojc+54Y6Lshbls1lf8sB/v8DG1TvSnsu2lXdf/oT3pi9ErSh46LMK10WY+cqnnH/NV7OWzytOx7Bno0hPEZkbtz7JtWIgIhOArao6T0ROi8vzF+AOHFP4HcA9wLVkYC5vTTxduarWAD0Stu3A8RZKzLsRGB+3/irwasvE9EnF4k9WM+2pf8dGMCVXBD7tAyO3nfVGCyoHK5ds5DfXPUaoNo0bcgJqa0bhOo08BQWMdQx7ZLuqjkqxbwxwrtvxWwZUiciTqnpFLIOIPARMc1c9m8tbE3/EcBGzftU27r31OaLxNauYV4tPxyW+vyCNk0B1bZi//L+pWJZN/yG9OOnMI+jdp2uTfIlsWreTe299rqkCSOeQEJPJMJxIsx4oyUN/QAwrB+MEVPVWXMcZtyXwM1W9ItZ36mb7JrDQXZ4KTBaRe3E6hocBH7dYkAzxlUARoqr87X+n8a8ps4nWht0WgNsR7NOxEWnaIW2aTc00ItRWR5j6VIPd/ZHfv8aVN32Ni649NWXxTzzwJs89/B6RVJP+JHNIiJcpttycIjAMqqtD6fPkiNiI4Vbk/0RkJE6jYzWuWV1VF4nIszjjpqLADfn2DAJfCRQln3y4nOnPfex09LnxbBxTgVsLiyVoXDNLtZyLfIm05rlaQ/ZCkylT2aHhGUgk5iQQf3ySfJGwxRN/fItRJx/CoEMObLJ/yfw1PP/oTCLxz12yc5lm+nPFOv9jiiA+b9x1dKqqIF/YOfIOiqGq7+IMoEVVv5Mm3504nkRthq8ECpSZbyxk8t9msHHtdixLnXdOFUSww9GGue7iY/r7ncA+qUj6MW764QvXRfjhBQ9glJjuYQblFUH6DuzJqiWbGrx6jDR9T16ew2YCCALUpgt7kkOcAHLtMy6QF3wlUIBMefg9Jv9tRgo3OgVbG/oH/c5gn2wwk7caFMBWrGismm+xb4/F0gXr6+d4cIKmtvJzZxrYeRrjogiRDhw2wlcCBcDe3TWsXrGV2pow4VCEp/7yjtPkToVhoJbdWBH4ncE+uSLVILfE2n8rP3edu+bHHKRKzgaLFSO+EmhDbNvmb/e+zrTn5hCNxfZXhWgzw1Lc+D6NzMLuX41bT7acbl8m+bycq5hl6oiy128wUjx9LXjuMpEptnHP3trkcuQczwPB2iW+EmhDXnluDq++OLdBAXgl1hRX6mtmiS+RplhOt6+t8hWiTB1OdqFhxrZUtOC58yxTTA6B8spS8oHitwR88kw4FOWd1z/nofvfIJJo9xepr22lbw2I+8L4ZiCfPJKn507Jb5Rzv2PYJ2/U1oa56dpHWb9mO9FQNPmH3u1w82O/+XRYDKEunDq+UC5R2u/8wV7wlUCeeeW5uWxYt6N+lK+SpMYf63CL9/fPFNedNGf58klHkD2ftPR+xm9LzJPrmkqsbBG6dMtTxzAQ8R47qN3Rca+8jZj59qL6CIuxATNJFQH4fv8+hUfi89gKz2dMrezZXZc2X+6QDj3RvK8E8kxFfGdX3BB63/Lj4+PiVn6Cpfn5PDn9D36fgE+eOPdbo1ny+XpCsQ5hDyMnfXw6IsGy/H2e/JaAT94Yc9pwJlx0PC8/M7tx9E+S+Er7+HRUBPbnK4CcStG2BESkEqhrSeA5XwnkGRHhBzedxQWXnsDn89ewfds+TNPgsYfeJVSTEJmxtUco5TJfkwstAJmKWfZClCkfsoMzVkGEzknmSm4NnI7h4miNi4iBM0Pj5cDxQAgoFZFtOPO2TFLV5ZmU6SuBNqLXAVWcMe7I+vXyylIevO91wnlyi/PxKVRiOqGmLkW46pyTuzmG88AM4C2ceQsWqqoNICLdgdOBu0TkJVV90muBXqeX7Ao8DIzA+Y2uBf4DGO5m6QrsVtWRSY5dDewDLCCaZlaeDs34846hrDzIo397hy2b9wK+ecinA2NK3gaLOR3DRdMnMFZVm4RXVdWdwAvACyISzKRAry2B+4HpqnqRO2F8hapeEtspIvcAe9Icf7qqbs9EMC/U1UX4cuUWunapoF/f7rkuPu+c8bURnPG1EQA8O2UWf3/kPcJhf6IYn45Jp6r8mIOgeEYMJ1MA2eSJp1klICJVwKnA1e4JwkA4br8AFwNnZHLilvLPVz7hrw/PwDQNopbN4IE9ufO3F9Kje6d8itFqnPfN45j53lJWLN9MJNJ4YFkqcy1p9nk1+Xopo1ACnhWDTMUse1vIFL+8Z19+AsjlesSwiJjAXGCDqk4Qkd8B38D5bq4ArlHV3SIyCFgCLHMPnaWq16cp9yycb+2DqvqpiFwXm+S+JXhRf0OAbcBjIjJfRB52e6RjnAJsSdMZocAbIjJPRK5robwAzP9sDX95eAZ1oSjVNWFCoSjLV2zhv37zQi6KLwhKS4Pc/6cr+fVtF3L2OUfRf0APxGj8Amvc3/gUvy8xHynyZVJGvvKl2taWMiWTpb3K3mbPgiGoKZRX5CeAHDgTzXtJHrkJ5+Me401ghKoeBXyBOw+xywpVHemmlArA5UfAz4ErROQMoIn5PRu8mIMCwLHAjao6W0TuByYC/+3uvwx4Os3xY1R1o4j0Bt4UkaWqOjMxk6sgrgMYMGBAWoGefXFOg5+9i2Upq9dsZ936nfTvV/ymIQDTNDhpzDBOGjOM/fvruOzbf0k572qylyrTfLkooyPI5DVfIcrkNV9by6SALely5g5ViNi5MQeJSD/g6zhTRt7slK9vxGWZBVyUZfHbVHU38DMRuQvHO6jFeLny9cB6VZ3trj+PoxQQkQBwAfBMqoNVdaP7dyvwEjA6Rb5JqjpKVUf16tUrrUA7d1Yn3R4IGOzeU5P22GKlU6cybr/9AqBpzcpPfmp3yYC6UL6ml3TGCXhJQE8RmRuXEq0b9wG3AKmmXLsWeC1ufbBrYXlPRE5pRtR/1cusOhH4R2ZXmpxmWwKqullE1onIcFVdBpwJLHZ3jwWWqur6ZMe6ZiNDVfe5y18Dbm+p0CccP4SVq7c1mX0ratkMPbg3qkok6owBFJGEGFeKutvd60MA0zQx0sVSLwBGHj2Ann26sHXrHudNSYfSYIjNRb58Uoiy51qmYpa9NYmdX4SqLvmbaD6DEcPbU3k4isgEYKuqzhOR05Ls/yUQBZ5yN20CBqjqDhE5DviniByhqnuTla+qLyes/9Gr0Onw6h10I/CU6xm0ErjG3X4pCaYgEekDPKyq44EDgJfcD24AmKyq01sq9EXfHMXUVz91av3ux9A0Da698mSmvrGAx6Z8SHXiwKtmME2Ds049jP+8biwV5SUtFbFVEBF+etPZTPzv57C17d9XH5/WIFa/2bs/PwHklJy5iI4BzhWR8UAZUCUiT6rqFSJyFTABOFPVqZaqaghnsBeu4lgBHILTqZwWERkF/BIYiPNtFacYPSpToT0pAVX9FGii/VT16iTbNgLj3eWVwNGZCtUcO3ZXs9+KYBsgNqiAUWow/f0lrNu4k1AWbpWWZfPWzCVs2rKHP/3PZbkWOWecMGoIxx8/hNkfr2y2MeDjU5SIkwLBfLlt5iZshKreitvp67YEfuYqgHHAL4Cvqmq9vVpEegE7VdUSkSHAMJxKtheewukk/pzUpidPFOWI4ckvfkzEstGgUf8hjNg2X67e1qJyo5bNwmUbmTJ1LhecM5KSYGHenssvPpFPP1+blbLz8SkWgqUZjXlqEa08x/CfgFIcxxhocAU9FbhdRKI4g2mvdwd9eWGbqk7NhXCF+ZVrhpVrt2O30nDCqGXz53+8x5MvzebBOy5lYL8erXKelnD0iP5c9e0xPPbkB0SiTSsBfgvBp+gR2F+bJ3OQQsTObewgVX0XeNddHpoizws4o3yz4Tci8jDwNq5JyS3zxUwLKkolcNjQA1mxehuW1aJWUEps22b33lp+9fupPHHfNc0f0AZc/q0TGX/Wkcz7bC1bt+4lEDAQQ1i4dAPvz1peP8CsKIKGFYJMvuxtI3uyEWUGIELnzvkKIFeU00teAxwKBGkwBynQMZTAt785mtffXdxECSQ+S5lS/2y6nkOr1u3gtMvvo1vXCi4adwwXn3MswUDrRhvcu7+Wx1+czWvvL6amNoxl2wgS6/YBEcT14igtCXLsEf257pIxDOnfE4Bzxx3Nt294hB27qhvdn0Jw/PDx8ULsPaypy4+LKLS6Oag1OFpVj2w+W/MUR8CMBErLgki5iUqDX7FtOL7FxG3LNAGN7ogC4UiULdv3MWnKh/z87pfQVpz9vaYuzDUTn+SZVz9h995awhELy1Kilk00ahO1lGjUJmLZRKI2+2tCzJzzJd//5WS+XOP0h5SWBpn0uys4fcxwykoD9fcjdj1+8lOhJwAMsFvxXYtHcbyDvKQCYpaIHJ6LgoqyJTBl2lzqwhHsksY6TKKKSuvUeCNRi0+XrOeZV+cxYnhfhg7oSVlJyzquVJXVG3eydec+SoImC5ZsYOee6owf/tpQhD9Pnsm9t14IQI9unfjNzRN47IV/8/cXZxGt68jzJvkUK52ryvJ2riKcVOZk4CoRWYXTJ9C6LqKFxrxF64haTT+UaoJYjmbPRZArNak3DQGEIhb3PfkeYkJpMMDNV53Buadl1yLbsHU3P/3dS6zbvAvL7eQ2bLJ29lq0fFOTbXMXrnU6jo3G5aYy37Zl0DAv+QpRpubyFaJMhX4/Y8u79+cpgJwK0eJTAuNyVVDRXTlAWSrXMRHUbNysTNbUjH8JkuaTpgogHlWoC0e59/F3+GzZhozlt23lx//zHKs37qxXALFys20AJ3Nn7dO7izMK2mi4L5Bw/UbDtozvU57zFaJMmmJfsv2FJnuqfW0uu+G8f/kctFmE5qCeqromPgEZtwKgSFsCtrihH5LtFIFA+oc7FZ72xZ20LhzlmemfcPTwvmnlBefD//6nK5g6cyEr121n8/amI8PVcFoymaLAzuoaLp74WKOaVDgcbTAtJbkv8edN3O71PrVVvmKQyWs+X/am+xSwPEvYMpSimlQmxkMicpWqfg4gIpfhTPT1SqYFFaUSqA1HnBps63iIpsRpHTTetn33/maPs23l5w+8zOyFawhHLLAV0SRKTMjqulQgqsrqTbua7BMTDI8zVubnlfPxaR41nPc8XxShErgIeF5ELsfpH7gSJzZbxhSlEujSuRwM5+MHkLpZQMu+bLFyYymBkqDJV0YOSVtEJGrx0YJVzFm01lEApJEVml5XKpniytEU8oGjuJLGGYovO3Z8srZ6svNmi9cycp0vFxSi7G0lU2vKHvdsd63KTwC5YhwnoKorReRS4J/AOuBrqppVJ0pRKoHdNXXOdyrmCZTu92ul39Y0hO5dKrlwbPLQSO9+8iW/n/wOW3buBytBDBEQddz+4zbHm5yava6EY1L2a4mghkKiImjmBfTxaStiz/Tu6vyMGIbiGScgIp/TuJrWHTCB2U7E5A7iHRQMmI5ppg1C55SVBOjZtZJzTj6ci8cdS+fKpm5sHy9eyy//9q+Gmj9NKz5quBsTTD/12zNowaTrxAbAFNTWvJvPfHyyQQUwIBhs3YGZ9edTiOZoUpk8MCHXBRalEigpNdN2dLYmv7z+bL52wqFp8/ztpY8aKYBEF03AbQ2krsHn/LoMSd1a8PEpQIIl+VECUFR9Amu1mRGrIiLN5YmnKJXAvtpQm3RimgGD7XuSz2oWz7qtCR20rs09TwMgfXzaBftq8jWfQFH1CcwQkReAl1V1bWyjO9fLycBVwAzg714LLMq6YefKMscKFpPea59ANvncpikGlJaYHNy3Z7PyDUmWxy2jSSdza8qebp/XfK0hk9d8hSh7IcrkNV8xyB5730yo6pSfAHIAquIpFQDjcAzhT4vIRhFZLCIrgeU4873/QVX/nkmBRdkSqIlEGzqG89RiDAYM+vbqwvGHDWg27w8vGMMNv3uOUML0l6m8jHx8fBqINZhrIn4AuURUtQ74M/BnEQkCPYFadwL6rPDUEhCRriLyvIgsFZElInKSiPxWRDaIyKduGp/i2HEiskxEvhSRidkKGo+t6nSG0nQ0Yrapkcw43j+9u3WiNBigsqyECWOOYNLESzzNQ3z00D7cf/MFDOvfs76/Nldy+slP7T0BqAlWvgLIaW5HDIuI6U4eP81d7y4ib4rIcvdvt7i8t7rfxmUicnZmcmtEVTe1RAGA95bA/cB0Vb3ItT1VAGfjND1+n+ogETGBB4GzgPXAHBGZqqqLUx3jhc6VpY7rYw7bMbddezZfPzEnQfkAOG54fybfdiUAsxev4eYHX6YusWXg4+OTkqrK0jydSbBy6x10E7AEqHLXJwJvq+pdbkV4IvALNwropcARQB/gLRE5RFXz+qFo9spFpApnGrRHAFQ1nIHmGQ18qaorVTUMTAHOy1bYGLtrahvXGjwsp9qnAAZs8TDyN1tGHzaAEw4f6Fmm1shHhvnyIVOmsheSTL7srZcvtrwrTx3DQM76BESkH/B14OG4zecBj7vLjwPnx22foqohVV0FfInzzcwrXtTfEGAb8JjbxHlYRCrdfT8WkQUi8mh8EyeOvjij2WKsd7e1iIqykvqO1nQfufjlVPswoKwsSO+unVoqVkpEhF9deRZm3MTZzX2gc50vU8WQD5kylb2QZPKarxBl8pqvrWSKvduVZfkJIKfk1Bx0H3ALjZ3CD1DVTQDu397u9oy/jyKyT0T2Jkn7RKRpQDIPeDGoBIBjgRtVdbaI3I/TnPkTcAfOPbwDuAe4NlHmJOVpkm2IyHXAdQADBqTvfLWImzcgBx3DQdPgzGOGtbygNHTvXMEpRw7mw0WrCUd9s5CPTzoUsCTpp6JVTpZB90NPEZkbtz5JVScBiMgEYKuqzhOR0zyU5fn7WL9TtbNnST3iRQmsB9ar6mx3/XlgoqpuiWUQkYeAaSmO7R+33g/YmOwk7o2cBDBq1Ki0NyIUtZxRhS18RgTo06OK+354HuWpwlPnkDuuHsdv/vEG7y1YQdSyWyq+j0+7Je8B5Lx7B21X1VEp9o0BznWdZMqAKhF5EtgiIgep6iYROQjY6ub3/H1sTZpVAqq6WUTWichwVV0GnAksjl2Um+2bwMIkh88BhonIYGADTifIt1sq9ClHDGL5xm3OhCmQXBkoqd0xlfqh6Z27lTHUg+9/LqgoK+F3101gxabtXHrXU0SsNPKnIt11FUK+fJ4rF/kKUSav+QpRJq+kKiM+gFyexglojjqGVfVW4FYAtyXwM1W9QkR+hzOI6y7378vuIVOBySJyL07H8DDg43TnEJF9NL17sXVV1aqkB6bBq3/NjcBTrmfQSpyZ7h8QkZGuAKuBH7hC9gEeVtXxqhoVkR8Dr+MYbh5V1UWZCpnIt087lpf+vYjte92pGDN8IOu/uQYs27CN2nCE8hZOFZkJBx/UkyvOOJbH3ukZCwkAACAASURBVHJblW3gohy7B8XhHe3TkYg9m7tr8zOzGGRkDsqGu4BnReS7wFrgW845dZGIPAssBqLADc15BrWVOQhV/RRIbAJ9J0XejcD4uPVXgVezFTAZXSrLePbWK3ji7XlM/2QZ2/ZWE45kMAWFNCQRIWDkf+D0T849mednfc7e/aG8nE+A0oDpmNLitvsmKZ+Cw303g4H8xQ7K9WhgVX0XeNdd3oFjQUmW707gzkzLFxEBLgcGq+odItIfOEhV07YkklGUI4YBulaWc+O5J3PjuSezc38N59z2CHURj7OnuBginDbi4Lw+bDFEhCtOO4aH35yTl47i0pIAd14xjl9OeZ3aGsfW6rcCfAqZkpL8fJ5Uc68E8sCfcTyQzsBxzNmPMybr+EwLKsrYQYl071TB3VeNp8wNPaseU0VZCb+5ZGybyAxw7ZmjOWn4gJy2RBIfZUOE8pIA917zDbburSYatVLOK+wnPxVS2lubx/kEim+O4RNU9QagDkBVdwFZ+dQWbUsgkdNGHMyM/3c9dzz3NtPnLyNqqbMj5kUkccs4sYCuOXMUVRVN5wPIF8GAyQPfP58Vm3ewYPUmtu+tRgRKAwEito1l2ZQEnDktw9EoAcMgYBqEohaoUloSIGrZRC2bA7t15qiBB7Fq605WbdlJwDCIWDaDenfjpEMHUl4SpOTLdQRMk7BlO2E3lAZv5nhvq3SeV/nMl2qf13yFKHshylRIssf+ClQlmaujtSjCCL8RNyKDAohIL7KcMaTdKAGAitISJl5wOh98sZq9NaGGSdbjcR+0irISLjrpyPwKmIKDD+zBwQf2yElZA3t347QRByfdN+rgfgzo1ZWlG7Y5G4ScjLPw8cklsbe2JpyZeTf78wl28UwqE+MB4CWgt4jciTPn8K+yKajorrw5ulSU8fRN3+arhw+mNGBSURKkV+eK+oaAaQhfPXwwT//HZXStzF+o2kJARHj0hxeBtH1T309+SpUAMMDO41R4GclWAKjqUzgjk/8X2AScr6rPZVNWu2oJxOjXowsPXJtZiKL1u/aws7qWQw7oSVnQuS2WbfPFlu3srqmjPBhgaO8edCrLV1Cr1qFzeRljDhvEB8tWt7UoPj5p6Vyep3dNKcaOYVR1KbC0peW0SyWQCTv21/DjyVNZsmkbQdPAspWfn30KB/fuwX9Omcbu2jpsVQQImibXnHwcN535FSTdnL4Fzq3nnc6E3z2GrU3NtSSsp1rOZ77YPq/5OqpM7UX22PKuuvx1DBdUNd8DIvI4cFMsmKcbu+0eVU0M3dMsHV4J3DD5ZRZt2ErUtgm5Jsj/m/4eihOeIoYCYcvi8Y8+YVCPrpx/zBFtIm8ueHH+QjQIGm68XeP+JvbfJVtOly/ZO5WvfIUoeyHK5DVfvmVXdwa+ytL8BJCDomwJHBUfzVlVd4nIMdkU1KGVwNqdu1m6aTtRu7HtsS6N335dJMqvp77FIx/NI2AI/bt35cJjRnDqsEEF0zpQVT5YsYbn5y1kzc5dWHbDK6fAqi07HWtrml/fa8UoXb5clNFW+XyZcpsv0zKU/AWQU8C2C+PdzQBDRLq5rqGISHey/J53aCWwY38N2Xy3w1Gb5dt2ALBky3beW76KMw89mHsuHN/mikBVmfjS67y26AvCVhJlpk4SwQnFnb++Nx8f7wjU5Wt6SQWKryVwD/CRiDyPcwUXk8XIY+jgSmBIr+5JRxkn2jGT7YsnFLWYsWwVs1et48Qhzc9B3BIs265XNKqKmTDQbN7aDby+OIUCgMYX5SoCp7BmThxvRM5kX6HnK0SZvOYrRJm8kqqMuLECXSvyOdF83k6VE1T1H25I6zNw7toF2c7Y2KGVwJoduwkGDMLRhkCyibZKSVgHkjrW1kUivL1sRaspgWkLl3L3m++xdX9No+0HVXXml2d/ldED+3PH9BlM+2xp8+9oXAtA4l46H59CoKFjOH8B5IqtYxjA/ei3aKpe6OBKoLwkgBkwwLYbm0XiasdNaghup1UihggVJa3TkfXm0i/5r6lvEEpSu9+0dx8/ffE1eld1YsOuPd7mWYi1AOyifPZ9OgIG7mj5fOBt6sj2SodWAkN79eCAzp1YvXN38pGzGTwXAdPkvKMPy5ls8dzzzgdJFUCMkGWxbtceZ8WLEojl80cL+xQwpXkKIAd06NpQh1YCIsKfv30eV/39OXbV1NV7CQUMg06lJVSHw41mAEv1nAjwo1NHM6Rn95zL+NGqtazasavZfPUmIHdodLHZOH18EtmTr3ECClok3kEicnO6/ap6b6ZldmglADCkZ3dm/Of3mb16HV9s2Y5hCMN792T0oP7sD4X4cMUaFm7cwmOz5mHZJB0BU1pictVJx+VcttpIhBtemJp5P5yB6/dG2wYDa2m+QpSpPcpeSDLF/gp0Kc9ncMfcKAERKQNmAqU439fnVfU3IvIMMNzN1hXYraojRWQQsARY5u6bparXpzlFbFKZ4Thho6e6699wz5sxHV4JAARMgzEHD2TMwQMbba8qL+OcEcM5Z8RwtlZXM33J8iZeN2WBAN89aVR9qIlc8uGqtYDUd+SmcxKJ/a3P45t7fIqU2PNcE81PALlGJ205IeAMVd0vIkHgAxF5TVUviWUQkXuAPXHHrFDVkZ7EVL3NLeMN4FhV3eeu/xZovdhBItIVeBgYgXO7rgUuwNE+YWAFcE38CLa4Y1cD+wALiKaZpLmg+Z9zv0b/bl14fPZ8qsPOUNtuFeX88JQT+M7xnn6/jAlHo4Si0fqO3JQ+/bGObL+j16e9YICVz0EsOXpxVFVxJngBCLqpvnR3RrCLcVw7W8IAnG9vjDAwKJuCvFZf7wemq+pF7jzDFcCbwK3uPMJ340yw/IsUx5+uqtuzEbBQCJomPzntK/zktK/k7ZylJcGG0cxeavZ+zd+nHVGVxwByuRws5sb5nwcMBR5U1dlxu08Btqjq8rhtg0VkPrAX+JWqvu/hNE8AH4vISzhX8E3g8WzkbVYJiEgVcCpwNYCqhnG0zhtx2WbhxLP2yQG7amv50dRXmLtuAyqKxD2g2Qboiu3zmq8l58p1vkKUPVOZiln2fMiUbHlnKH8B5DJwpOjpDtKKMUlVJzUuSy1gpGtBeUlERqjqQnf3ZcDTcdk3AQNUdYeIHAf8U0SOUNW96eXVO0VkOnCyu+kaVZ3v+Sri8NISGAJsAx4TkaNxNNxNqlodl+da4JlU8gJviIgCf0u8YT5NuWHqK8zbsAEbdWr3CuJ2RSQOZkvclrgv23yJ/XepltPt85KvEGVvDZm85itE2b3IlEvZ6wPItdK4m6R49w7a7tWkraq7ReRdYBywUEQCOGb04+LyhHD6EVDVeSKyAjgEmNu0xCblz8P5HrcIL0ogABwL3Kiqs0XkfmAi8N8AIvJLIAo8leL4Maq6UUR6A2+KyFJVbdKLLSLXAdcBDBjQuqEXCpVt1dVMmjuHj9evb/ySiEIg/UPqtSKT63y5oBBlbyuZfNljfxXN46QyuYpV507zGHEVQDkwFrjb3T0WWKqq6xPy71RVS0SGAMOAlWnK/0BVTxaRfTS+5QKoqlZlKrMXJbAeWB9n13oeRwkgIlcBE4Az3Q6RJqjqRvfvVtd+NZokrkxuC2ESwKhRo/L5DSoIlmzbxsXPTaEmHGloqrsdwk6Hb4e7JT4dGYEaK58B5HJW2kHA426/gAE8q6rT3H2X0tgUBI6p/XYRieI4z1yvqjtTiuooAAGOUNW1uRC4WSWgqptFZJ2IDFfVZcCZwGIRGYfTEfxVVa1JdqyIVAKGqu5zl78G3J4LwdsbE996g+pkURNjHcIal7y0XPOZL9HA7CVfa8vkNV+hy9TRZI8bK9CtPF8B5IRcdQyr6gIgaVx/Vb06ybYXgBcyPIe6FeqcDE7y6h10I/CU6xm0ErgGmIMzIOJNN6rlLFW9XkT6AA+r6njgAJyOkdi5Jqvq9FwI3p6ojURYvG2rs+IOlFFNeEfc7T4+7Z1YpXxXHjuGi7ChPUtEjlfVOS0tyJMSUNVPgcTOkKEp8m4ExrvLK4GjWyJgR8A0DCc8dMyiZuKYgPxY/z4dEdcMGgwkCdfbWhTfu3Y68AMRWQNU09AncFSmBfkjhtsYVWXR9q0M6d6NL7bvaNhhkDRktY9PR6E0kKfPk5Izc1AeOSdXBflKoA2pi0a4atoLfLJpIxHbRnx7j49PPXtC+ZtPIE8zWeYMVV2Tq7J8JdCGPDD338zfvJGIGxhIDQWl0eAwH58Oh9v/1TVvHcMUXZ9Aimiie4B5rvneM74SaEOeX7qIcPwk9zF3UFG/VeDTYYm5Q9dEw83k7NCMctMr7vrXcZx1rheR51T1/7wW5CuBNiSSOFGMPy7Ax8fBgGgePSOKzRwE9MCJIrofQER+gzOG61ScUcS+EigGzh4yjGeXfN7YMcEPAe3jA0CXsjzNJ6BkEjaiUEiMIhoBBqpqrYiEMinIVwJtyM9PPIVXV37B3lAdhRc2LB/hxQpRJl/2tpE9/lzO9l1hf6L5NEzGGSvwsrv+DeBpd1BuRpPP+06IbUiP8gpu/cqprj90spBaseV0+7zm01bK195k8mVvG9njzmUAJlSWBMkXot5SoaCqdwDfB3bjdAhfr6q3q2q1ql6eSVl+S6CNOXfoYfzv7PeIkFELzsennaP5dd0voA+8V/IZRdSnFakMlvDMNy7lh2++zOq9u/E7hH18cALI5dM7qAO/dr45qAA4rEdvZlzyPX4++hRKg26vsKRI6fa1VT6yyJd3mRSMWLu+2GQvAJnyJbvRkLqV5WecgFdTUCGZg3KJ3xIoEESE7444jhdXLGLdvj2ErNgk2/Edcz7Z4b69rvttdvc03THxXwf/t2o5zv3cHcljx3CReQeJSClwITCIuO+4qmYcpdlvCRQQZYEgL3/jCg7p3p3kHXN+yiq5NUyJ1TQlm3tKo2VJWG/49hfA9RZ9cn6noJG/z1MRtgReBs7DmdCrOi5lTIdqCagqn+7YyJJdWxjYuRsnHTAIQ9q2BqCqLNixiUW7NtO/UzcO69qb5fu2IQF134fiqqEUBvFva+P7J9JkU0YEDZPxg4ZzTK8+3DVnBnW21fxBPllRnq8ActD4kWkBIlKGM2lWKc739XlV/Y2I/BbHm2ebm/W/VPVV95hbge/iTCrzE1V93cOp+qnquFzI3GGUQG00wlUznmbxri2oKoYYHFDeiSlnfYeeZZVtIlPIinLNjCl8tmMTimKKQedAKQHDJCSW+2AWVvWjaBByfP+EMjPAsK49uf2Esyg1Tf668N9srt6fo/J9EtkTydN8Armt5YeAM1R1v4gEgQ9E5DV33x9U9ffxmUXkcJwZx44A+gBvicgh7mT16fhIRI5U1c9bKnCHUQJ/WDCTz3dsIhRXc1u7fze3zn6Vh776rTaR6Y+ff8D8HRvj7P9QF40StXA+YvEzigmNWsuNkLjtkiJPJvsKPV+qfYmdjrH8dpJticsx4u+10EiZPDL2Qr5y4EBEhEeWzGaPVY0EbdQSsBMFyFD2tsxXSDLF/gp0Kc3TiGHIWV3BnWY3VjMIuild6ecBU9wJ51eJyJc4U/D+O1lmEfncLS8AXCMiK3EUj7in9+cTSMULqxY0UgDgxCZ5b+MKQlaUUrP1bsWCnRu5b+G7fBr74KvjCFoXafpsWNiOzSLmJJ34YfNJi/ONce9rzNTX0jAcoogIG2t2c/XMD1m3fzfba2oJ2RYiIAFFFXcSoJgy8GkZSq2Vv7EzksMwRe78wvNwJt56UFVni8g5wI9F5EpgLvBTVd0F9AVmxR2+3t2Wigm5k9TB05dPRLoCDwMjcJ7ya4FlwDM4vdOrgYvdi0o8dhxwP86r+LCq3pULwTMlaif/lRVFNUfVgCR8uGUl173/LCE7mmRvqpqjEhQhYvumoEwQEQzDxlZF7cRR2C0pWEFsbp07rf7nsqONfzun01lRO11V18crYihWYd7HniIyN259kqpOis/gmnJGut/Nl0RkBPAX4A6ch+MO4B6c72iyGkPKC8/lPAIxvFZ/7wemq+pF7jzDFcB/AW+r6l0iMhGYiDPxfD2uRnwQOAtHw80RkamqmlFsi1wwtu8wXlmzuFFkQgGO6tGHskDrDU+/7ZPXUygAYg04Eu0RImCYSqkaRPyOR88YgG3YiC2ImfsPSKNfKslvJwJiKkExMYQmLU+fzOhWmOag7aqaONVu8iJVd4vIu8C4+L4AEXkImOaurgf6xx3WD9jYXNki8jhwk6rudte7Afeo6rWeriKOZn2wRKQKJzzpIwCqGnZPfB7wuJvtceD8JIePBr5U1ZWqGgamuMflnYnHnEHP8koq3A9+mRmgc7CUu0/4equdM2JbrNq3I+V+MWJPXuO/YipRLLqVltXL27hDILFzIHE5n/mS7fOaL1cyQYlhooaNYeD28mmSfNmfS0wl3pGs4bdrTKkZ4NxBh3N0jz5xv10irXXfMy0j1/ly9Sw423dFsvJ4zBzNnYuoiPRyWwCISDkwFlgqIgfFZfsmsNBdngpcKiKlIjIYGAZ87EHqo2IKAMC1whzj5XIT8dISGILj1vSYiByNY+u6CThAVTe5AmwSkd5Jju0LrItbXw+ckI2gLaVXeSfennA9U1cv4rOdGxlW1ZMLhhxJl5LWG5UYEIOKQAnVKYa/iwCmgrqVShRx/dnLzSA3jhhDt9JKbpn9CjWRcJJOtVQva6p9yZ7idC98S86V63ypZT+mx0Fsj+5jQ/UeQDBMUFWnw7alMknDbxJP7Lc7qLwzR3Xry866GgZ27saFQ47i+F79sVWZsfFL3t24gq4l5XQtLWPe9vV8tn0jm2v3N3/eFu0r7mdBDAVxQqrkDQ8feI8cBDzuWkEM4FlVnSYiT4jISPdMq4EfAKjqIhF5FifyZxS4wYNnEIAhIt1iJngR6U6WfbxeDgoAxwI3uh0c9+OYfrzg2d4lItcB1wEMGDDAY/GZUR4IcsnQkVzCyFYpPxER4TtDR/HYF7NTmgZifuuScFsChsG5g46kc7CUHaH93Pnpm64pq6VPa+6e9vyTXPZDe/bkX+u2NjLRxDpsc3/+hke61DSZOPJMJgw4oklOU4Sx/Q5hbL9D6rd9D9hYvYezp/+VWiuSY9myoZCfBW3yTrTy6XJTjOoCktTIVfU7aY65E7gzw1Pdg+Mm+ry7/q0sygC8jRheD6xX1dnu+vM4SmFLrInj/t2a4lhP9i5VnaSqo1R1VK9evbzKX/D8x4jTuHDwSIISf6sTzRWJCf4y5lt0DpYCcMXQURzdo4+H4xKb117OVdxJDGVHXTWGSNxAsGzLS0X8/oblMQcOSqoA0tGnsguTTrmEnqWVGFl7EbX9fW/1JFBj5yeAnOB4B3lJhYKq/gMnbMQWnG/vBar6RDZlNdsSUNXNIrJORIar6jLgTJymy2LgKuAu9+/LSQ6fAwxzbV0bcAZFfDsbQYuVgGFw+3Hn8PMjT2dDzR4eWz6LqWs+J5LokaQ4LQKBg8o7cUKvgfW7RISfHHEqP/zwGWrixhTEjomnxBDOH3gUVw49gTc3LuOvS94nlMIzKlUZSfe1Vb4UxMwzlcESzu57GB9sXenU/gXHXTP+9no81wFllfz1K5cgIoh7gjnb1nDvohnUWpH6sgHKTZOfHPHV9EKm4KTeg/jo3P9g9b6d1ERDGGKgKO4/TAQLRVECYmCpuq1F4b8/mcaCnZuaXlcL72dGZbRiPon7272kwoPQOUC92fsLkE04/QdlOF5Lp6rqzEwL8WpDuhF4yvUMWglcg2vvEpHvAmtxmiOISB8cV9DxqhoVkR8Dr+O4iD6qqosyFbI90LmkjENLyvjFUWN5f8sKdoSqk9Y9y8wAtx379fqPUIyTeg/m6J59+WTHOiIpPuoGQreySn525Fi6l1bQt7ILT6yYTSRcS3v0XQ8aJsOrenPuwCNZsmczT3w5BwutVwaZUGYGuHv0uRzZo0+j7Yd1PYDXNy5h0e5N1FnR+v6aMQcM4chufVKU1jyGCEOqemR83J3HTeDSGY9Rm8rjrF3gvBl7IjX5PmXRICLfw+mb7Qd8CpyIM8DsjEzL8qQEVPVTnJntEzkzSd6NwPi49VeBVzMVrL3SvbSSV8b+gIeWfcS/NixifziEopSZQY7q3ocbDjuVo7s3HStiiPDIyZczZeU8nlwxh+21+7Gd+iOGCFXBMsb3O4LvDz+J7qVODWpvpJao1lF0T3izCH0ruvCdoaO5fMgoTDFYVb2ZYMDCaShlogGEE3sN4pYjxzKi20FN9hoiPH7KFUxeOY9/rlmAaRhcPGgkFw3OyhGjxRzW9UBeGvt9/rLkfT7YspJaK4Kq2wKKZVJncFvYsohSnG6qhigleQwgV4SvyE3A8cAsVT1dRA4FbsumoA4zYriQ6FFWycSjz2Li0WdldFyJYXLl0NFcOXS0p/xPrJyFpRamAd5tAoWKIqKYBgTFZFz/YVw77EQAVu7bxvyda1GxKQmAZacYhJeAAIM79+DxU65o0vKKp8QMcPWwE7h6WJs4tjVhSOee/G70N5vN99Cymdy7eEYeJGoNNK8B5IrQHFSnqnWu6bJUVZeKyPBsCvKVQDtm6e5N2O6H0yrykawCGO53OqIWS/Zsqt+3av92goZJyHbMNYYoqQdbO4UExKBbSQUPnnhpWgVQzOyMVCNio3mdpzE3CLAv6puD0rDeHY/wT+BNEdmFh0FmyfCVQDsmNgZCBEzDMR2pnfx5bxTjK64DtEm+uH25zpdMJsT5+Md/p4NickSXBnv8kM69Go2sNgwQVWJdJ/HX1aOkgquGjuHQLgdwYq8hBPJpcsgzw6sOpFMgSI0Vce6F19+EtnsWoOH3bs0xPI3QwvL88YKqxpqCvxWRGUAXYHo2ZflKoB0Tthv6A0TE+ai2NJhaXon/WjRogaBh8p2DT6pfH9ypJ8f3HMzsbSuIuGFBRMBMuNZSI8Cdx53PqQccQkfg7D5H8Melb1NXF0GK6ncHUEKavwByxdYSEKf5ejkwRFVvF5EBwEi8jTZuRPutBvkgBgTMWBWnAHy/M0yGKAHTwhQ7bptNl1KhMtB4NOl9x1/CpYNPoMxIHqqhf0U3fjfqWx1GAQCUB0qYfMp1nHHgoUnGJLT975suBUwbW/NXPS/CmcX+DJwEXOau78OJ05YxfkugHXNOn6OZv2MNISlWd0LnrQsGbRomBYAau44X1s7m6oNPq99Waga5ZcQ4bhmRk8mW2g0HlFdx3+jLGm37w5LXeHrVv4lS2DaQLiUFGUCuUDhBVY8VkfkAqrrLdeHPGL8l0I4Z1+doepRV0vCExz/pmrA9cV9L8uXqXGAmqX6F7Sizti9vst3HG98behp9K7sTqB/FXmjPgrN9bx4DyHlOhUPEjU+k4ASug+y0uq8E2jEr9m9ib3Q3hsTHHGqNj3XuyxAUUxpH7owhCAdVdG26w8cTnYPlPHPyj/n1Ud8k0OgGF8azYIjz21cGSr1eUosQitIc9ADwEtBbRO4EPgD+J5uCfHNQO+a5tR9iEcUQw1UE6VwFFbfruMAqPE0pMQJ0Kynlp/P/Qp0V4eBOfZjQ90QO6dyvrUUrGkrMABP6HsOm2p38Y8V7hLTQTIaK5NFlp8A+8M2iqk+JyDycAbsCnK+qS7Ipy1cC7ZgtdU648YBhY9lufJoUCEpAlOFVA1m8Zz12s+3fZB2NrY1QGSjFlBAvb5jpygiL967h9c1zuHbIOVwy4LQ8yNF+uPbg09kdruHFdbOxchKlNjcYooRs3zsoHaq6FFja0nJ8JdCOObHHcBbuXkMEi4BpNw2s5hKLtXN4VX8eOuGHbKnbw87QPv5r/hNsDjmKJBaaIGY9CBomtx95Bb3KOnPnwqdYXb0dWw1UJWXcnlgZicvp8tXLCHyl12FEtZrP9qxsckzYjvLoytcYe8Cx9Cit8nB3fABMMfj54d/g+wefwYXv30WN1TRyZ8y3v8HHX+sH2MX7/8evx5YzfRbi/3bLVwA5KEolkCv8PoF2zPn9TqB7aed690ARZyBVYhKBMiPITcO/AcABZV04rEs/fnnktygzAxiG43Nfn9cMML7PcXz1gMMZXNmbLeEdBE1182nSc8SXkbicLl8slZgmm0Jr+WzPipTXa2AwZ+eyvNzb9kbX0kp+dvj5lLu/d+Lz4fxVDEMxDSeER/zz5FQQtNnfsblnIVYOKPuj+/Nz8R77AzzOLFYmIh+LyGciskhEbnO3/05ElorIAhF5KW72sUEiUisin7rpr617sU3xlUA7plOwnL+feBNXDDqN3qVdqDTL6FVSxYk9hjOgvCedzDK6BSs5vfeRPHTCjxnRdWCj40f1GMpfR/+Ik3seRpdgBZ0CZQztdBC/OPxCbjnsAsCpSao6H4SgYcV1QucuBQRKzBDbQjvTXm/YjlBq+I3bbDmnz3H8/thrOaRTLJBe499BUAKGRdC0CRgxt93Y/tz97oISNCyCRh5HuHkVr3lCwBmqejTO4K1xInIi8CYwQlWPAr4Abo07ZoWqjnTT9ekKF5En3L83eb+49PhvTDunKljBD4aN4wfDsvOfH17Vl7uPuTrlftNwJlW33ciVAVHIceRKJ3yE5dQ8UdcPrqmdwUbpEsyjCaEdcmz3g3n4hBu58MPb2JPGRdMQpcRMdDZQkv0u2aGUJQ75bkVy1QetqgrEmjBBN6mqvhGXbRZwUZanOE5EBgLXisg/SLjhqpq+ppQEvyXg0yLe3DwLqR+MltsWQCyZhlVvKzbrJ3hvmq/CNNjSTGvBp3lMw+Suo79HpVnmeoylqxanW88+GSg1Vv4CyOXSRVRETBH5FGfGrzfjZmWMcS3wWtz6YBGZLyLvicgpzRT/V5wYQYfizPcen+Z6k7AxvhLwyZodod1MWvEcIkrQcEwEhtucN8RZNhKWJbYtST4h/hjHVzxoWI0GjDlmp4ZzOflsu/RHywAAIABJREFUd5swoOLANrwj7YdDqwbw4im3cXH/UzBwfgtTbDcpAcOmxHTue2w96CYjzW/c3LNgxp4lU6kK5i+AXAbmoJ4iMjcuXdekOFVLVUfiTPgyWkRGxPaJyC9xJpR/yt20CRigqscANwOTRSSlZ4OqPqCqh+FM0DVEVQfHpSHZXL5vDvLJmg+2z0fVeYltd/pDw4yvLqVaTiT1MSZOJ2RUDYjr4BZodC5TDAZUHsjhVYOyvRyfBEqMAP0ruxM0nakuGywPjZVyg5J2/gaSWnEyfRaUSGG6iG5X1WQTbDUtUnW3iLwLjAMWishVwATgTNdshKqGcPoRUNV5IrICOIRmavWq+kMRORqItRxmupPcZ4wnJSAiq3ECFFlAVFVHicgzQGwSg67Ablf7NXtsNoL6FB5hK0IUC1OcpnJuegIaAhmbaL35J4DtKoLkjO5+GLcc+p12OzdAW2GpRZlAnaYbZeIgNJ/HCwKYYmPnaVa02IjhnJTlhG+IuAqgHBgL3C0i44BfAF9V1ZqE/DtV1RKRIcAwnCl8mzvPT4DrgBfdTU+JyCRV/WOmMmfSEjhdVbfHVlT1kjiB7gH2eD3Wp33QrbQT0FAbNFv4CQiIiSEWiIVlC5Y2zBBmCJQk7b1zTBF19i4qAnkMONZBGNX9CB5f/TJB27t7TNN82XUYVwXyZA4CJPUsRJlyEPC4G9fHAJ5V1Wki8iVQijMBDDjTQl4PnArcLk7HmgVc77Fz93s4QeSqAUTkbpw5hltVCSTFjWt9MVlMcOxTvKgq72/7CAMbu9ELHm8ySDQfSMJy46NsjTj+58RmB5M4M4QmHO/g2JNhZfW63FyYTyP6lPfi/L5n8Oy61+O2Jvsdk/0+8euJ+9OXIUC1lb9xAjlpwgCuSabJBNSqOjRF/heAF7I4ldC48W2RpWuW145hBd4QkXlJOkJOAbaoaqqwjumO9SlS3t76Hsv3ryTg+oxL/ZuU6DseP59B42UDOLCsB+WmUBYIU2JG3YFCTusiYMTKjjmFOh2TsQ7FgGG7A5ec2cZ8WocrBk6g0jQTfmNo/Hs3/D6x399ZVwJuh6+keRbij4n97hVmfgLIQVEGkHsMmC0ivxWR3+K4nT6STUFeWwJjVHWjiPTGac4sVdWZ7r7LgKezPLYeV0FcBzBgwIAMLsEnn0TtKPN2zWPy2meBCBDAENJ0CCeuNyx3CpTyvSETmLTycepsJ6xFyG6olzTuAI7rBE5Svin5sR93RESE8X1O49VNM4hovDswSZaT/T5OHjPtcW38mxbWB75ZVPVet9P5ZJzX5BpVnZ9NWZ6UgKpudP9uFZGXgNHATBEJABcAx2V6bJJ8k4BJAKNGjSqyn6RjELJC/O/S/2V9zXpCtvNxKDGihG2vdQnXvo9B77Ie3DL8OlZUr6ifQUoESgyLsJ1Zrd7ERolkdIxPZnx7wLlUR2uYsfXfWLR+oLmAWIS1rlXPEU+B1fI9oaqfAJ+0tJxm314RqQQMVd3nLn8NuN3dPRZYqqrrszjWp8iYsXUGG2s2YmEhGChgiFBqRJ3GvNIQaI4k1l+FXx3xc7qXdOWgsl6ICGbMNcPNbIjWl5d4bOwchjSUHZuUfEjl4Hzcgg5LwDD54dDLuXzgefxiwa3sCdeirltw7PeJ77qJfxaMWKA497/Y75gYQM4ZGe566whUBSrzd4FFqARyhZcq3AHAS26PdgCYrKqxWe0vJcEUJCJ9gIdVdXwzx/oUEbbavL5lOhG3xm2KTVSdiY3qJ7FP2y2llIpBt5JK+pT3rt/ar6Ivx3c7lg93zKovoP6D4NJc2SVGyf9v78zj46iufP/9VfUiybLkTZY3jA3G4AECwWbJECAmIQN+4QUyrMOQISGTkAyPzMtjhoQweVmGfELIAglkApmXDAESQpKBIQ8IJLyJ+cCwLzab7QBh8YIty9ZiWe6tzvujqq2W1JJacqslte7Xn/qo+1bd6nNc1X3q3nsWLlx03kjUcgyTzmw7RoaEnyOTv/4UXDABWJGLZlFTKWuX4S9yd66zbHIP9XEVLF0w7hjSCJjZ68ARA+y7qEjbZmDVUH0dE4s7376NjkyPF7AniJMja15J/uNhJKhPwutfBvWSAz/O0zufJBUM18FBLJ92BOfsdyYL6uYPo59jpCS8BIEF4WL8gNe/d4uiBV8hcvS4/Q5EWNsiIFahBHLljBOoJJKmE8YV7PWNLrbeOhQuYtgxILuzXbzc/hxvd7/Bwy0PRW6bvSN34yU+Qgkxu2Y2Tcmmfvs8eZzY9Oc80vIIuRIDhDw8ltQfyOcOvrRkfRz7TlOyieaaZjZ1bwJZidd/sEX9gfvU+hX8eSpWaGMcI+kTwGcJU1M8DxxHGCcwbFd9ZwQcRVmz80l++sYN5MiSMyPAx4tWAoY3chYeHg3xBi5dMvAP9rn7ncvbu9/mza43I0NQ7Esp4orjy2dqfCqfOvBTw1PKURYuXXIp16y7hl3ZXWQsw1AT6orSf2SHkarMw+iuVD0BJuRI4LPA0YRBZyslHQJ8ZSQnckbA0Y/OTDu3vnkj2Wj+f+90b/TkHxhRINfQxORz3sILOLHpJPxBfPlr/VquWnYVGzrX890N15CxHlOTXyhcXHcAJzS9nxmJGSxrWIYnl/9wLGiuaebaI67l4ZbV/OKt28hZFqRwMbjguHzQX35h2LNcFAU+sNkIU0aEsR/1lVoYLmOwWAXZY2Z7JCEpaWbrJB08dLf+OCPgAGBXtp0/bL2LNe2P0pHpIlMwPJbAs548/p7CL/dQ+PJZUr+UlbNLG6FK4uCGQ1g197/x+233kw56Sh3GleC8hRdwQP1Bw9TMMRr48lk5+2TWtD3Jus6X6FkMztM/MlwiqjdRCkbOussjbAlMwIXhjVF1srsJ4692AptHciJnBBx0ZTu4fsPldGbbCL98HmEwec8XOKaAnIlciUN6D58TZq3kIwuG77Vz+ry/pNav48Gt97Ir28ncmvmcs99fOwMwDvnMkv/J5c//LWnr/SvaP51I6eQXhiuZC3CiGQEzOzN6+WVJ/wk0EtYZGDbOCDh4dPt9dGU7yD+9eQqgT8bO/FNcrISF27gSfO7gf2Ze7X4jkkcSp8xZxSlzVo2ov6NyJLwkMxJ17My00d/rZ+RzLHElWD79z/dJtpLJBy5MIKKcbRcAB5jZVyUtJCxn+eRwz+UmVSc5qVw3j7Xe1yttr6fwyb9nsrQwCdjAof5CJJTkvU2njNgAOCYeOcsv4Pa9T6xPe999xe+thBLMr9uf98xcOYpS92YC5g76AfAewrQ9EKbrv3EkJ3IjgUnOXZuuI5XL15LteZKLewG+GZkg7+FtBbO8Pa8bYg28b/Z/Z0dqG5JYPv14FtcvraQKjjFmaryBdLCNrClyGFBUKYzQiaDg/lEUNBbzDJmRwwsjhgmri0liZrKe/7HkS8S8SrqIVu6jysSxZnaUpOcAzGynpP5BOCXgjMAkYUd6M8/tuI8tezaQCVJ48mnwZ7Oh4yliCkhbf88dT0bSH/jbEVOcCxd9hqVT+9USGpCubBtr2x6kNbWR/eoO5c8aTyLuuToAE5kTm87g3s23IEsx8K9p+OMfJ0uaGPn8EF5BHqKkskgQWAdbU39ifm1l1oAmaLBYJqpZYEC+OM2IVjacEZgErGt/hHs2fYsc6V7tG20dAfFw+sdyZBlOhKY4ovG4YRmAd7pf5fY3ryBnOXKWZl3HIzy6/edctPh66mKNw/hsx3ji2Bmn0JZu4dHt9+LJJ2dZFtYtZVP366SCsIiWMBLKhUkCLdfroUNAPDIA4XuPXZmdUKmaMmblLCpTKb4H3AXMlnQ1cBZw1UhO5IxAlZMN0ty75bv9DADkJ3/CQXrMM3zL9k7yRc8EUeHrMF4gzsnN5wxLlt9s/jbpoMftL2N76MxkebjlNk6d+3fDU8wxbpDEqXMvYOXsj9Ca3kpjfCa1fh3Xrz+frtweMKIHjHyVOCNJ7ySBhZ5AWcswv67CU4plsgGSagizJCcJf19/ZWb/W9IM4BfAIuAN4Bwz2xn1+QJwMWFhmMvM7IEipy78DEWf8QzwfsL/wjPM7JWRyOyMQJViZjy74995tOVW0kGKYvlaJPAtRy76ghZmchyMmBIc3ngiM5NzS5Znd7aDnelN/doDsmzoeNQZgSog6dcyr3YREI76AnJhEXpBEARRGaHCJIH9YwnAaPRr6cpsoz42rWKyl3E6KAWcbGa7JMWBRyTdT5hy/yEz+4akzwOfB66Q9GeEiTgPBeYBv5e01MwGdMMzM5N0t5ktB9btq8DOO6hKeXz7rTy67Sdkgt2DHhfzjLiy9K4INvBWH2vkQ/M+zenzh/ej7Ss2oBeer/iwzuUY/8S8BFYQOxBTQEw5BrvPRHhMNtjKL9+8nLYiDw2jghGtYJewDXWqkLy7VDzaDPgwcEvUfgtwRvT6w8AdZpYysz8BrxLWXBmKxyUdXbqSA+OMQJVhZmzevZYnW+8gR7qkcJ18Dv+htjrf4xMHfJMjp5887JQNSb+OhXWHoT63XEwJjpx+6rDO5Rj/zEzsR31sJoVP/n6/+ywXbeH7hJclEa0NZC3N062/qpzAQz//5KeMZkl6umDrVzJXki/peWAb8DszewJoNrMtANHffD71+UBhgeyNUdtQrAQek/SapLWSXpC0drhqg5sOqioyQTf3vP0PbO1eHz20qGchbpBFX58A0BBhYGLl7IuYlpgzYvlOn385t73xD3Tl2jALMGBh3eEcN+usEZ/TMT6RxFkLv8Ttb3ye7lx7n7SDxZ+ohRGL9hkB2/a8WgFJo88ufTpou5mtGOyAaCrnyCitw12SDhvso4udogQ5TivhmJJwRqCK+K9tN9GyZz1GlsJL6wmSliNfsj1c5NXe7I6KCoH4FpZ/t173YHhcjZfkqBkf2if56uMz+NSSH/Fm11raM1uZU7uE5poD9+mcjvHLrORCPtD8Vzz0zg2kgp4VAQNyBSsCYZKSnkRzEHoIza7gvTEa3kFm1hbVAT4V2CpprpltkTSXcJQA4ZN/YWTlAkrIAWRmbxarJwC8OVw5nRGoAnZltrKm9VZear8XyC/4BlGen4LhOFYwHuh/0/vymR5vpivbStZSe4+LqYZjZ52Pr32/XSSPRfWlu5U6JjYvt/0GnxRx+VGkQFiFbqhSRL7irJh5dkVkLGcW0chfPxMZgFrCErzXAPcAfwN8I/r7H1GXe4CfSfoO4cLwQZSQ+qHi9QQkvUEYlpwDsma2QtKXgb8FWqLDrjSz+4r0PRW4HvAJy05+Y7hCOgZmV+Yd7n7zY6SDXVDghucrAKOEhG/54u8xDmk8mfc1f5rN3S+xeutN7ExvpM5v5OhZ5/Pu6WcMcR6Hoz9ZS/eqQhaUsAzpYTT6UOPXV0DCfLBY2UYCc4FbokAuD7jTzP6vpMeAOyVdDLwFnA1gZi9JuhN4GcgCfzeYZ1ABY1JPYKWZbe/T9l0z+9ZAHaL/iBuBUwiHPU9JusfMXh6+qI5M0MWu9CZylt77PPX8jlsiA2AFBV8UJXwLiA0RRNgQn8cFi29DBY7ai+uPYXF9KQ4KDsfgHNRwMs+03kaOdFSFbLD70fAJopoVOV7e+SuOmnVxZQQtUxZRM1sLvLtIeyuhT3+xPlcDVw/zoyZMPYFjgFejWsNIuoPQJcoZgWFgZqxtvZENbT8jiAq95NkdJMg/zcfJkSJGT2hXMT/sEF8JPHw+MPfKXgbA4SgnR0w/i1c7/pOOzBaytofe92T/1/mHloAMm3Y/wVFUxgiUcSRQKSpeT8CAByUZcJOZ3Ry1Xyrpo8DTwP/KR8AVUMz96diRCDqZ+WP7nWxou6OfAYB8Mreeef+kZcnuzciivccIj/l1RzOv7gja0xtpSMxjWeMqpsRmVlIVxyQj7tVy9v4/5LVdq3n4nWsjQ9ATGQDRojBBgZNCSEzJyghZxjWBSjEW9QSON7PNkmYTWp11wL8AXyP87/sa8G3g4336lez+FPnbfhJg4cKFJYo1OXhl508ISBXdF1eOlIULbhCldCgytk169Zw2/2vEvAp9sRyOCN9LsLThFNK5HTy7/Ud7DcHgGHGVMjVeDiZk7qC9mNnqfelfUsSPmW2O/m4jTFp0jJltNbOchWGBP6J4lFvJ7k9mdrOZrTCzFU1NTcPRoar5U8ev2JNrHXB/TAFxsgzse+3TEF/Aqv2+7wyAY0w5dNrZHNJ4Br4SxFRH33QRhVuCLGaDR7uXFbPStipkyJGApCmAZ2ad0esPAl/N+7xGh50JvFik+1PAQZIWA5sIc2T8VXlEr37SuTZebL0WEWCD2OuEFxC3NAHiuOarmRKfg2c+gXLU+I00xPdz8/6OMUfyOGb2pRw58yI6M5uZEm8mne3kvrfOI7BMOLEp4WH4SjCn7rjKCGYTr7xkOSllOqiZMOotf/zPzOy3km6VdCSh6X4D+BSApHmErqCrzCwr6VLgAUL/xR+b2UujoEdV0tL9BMInrhRp25vDs+ixklGDR10sSXPt4ZUT0uEYJgm/npl+mCW0xm/k8Jkf45Wd/0YuWjj2iJPwp7J02vDrU4+YKn3KL4UhjUDk2XNEkfYLBzh+M7Cq4P19QL/4AcfQeEoghCdIkCVjftFJn3xh7oQS+JVaTHM4ysShMy6mMXEg69tuJ5Xbydy641k2/aMk/cplEZ1oC8PlxEUMjxOyuQ52dD8MwPTaE4j7jcyufc/e/Z4gOcRCmefFmVlz1KjK6XCUg0yujdbdq+nO/omE38S0muW8f8GPxkweBZN3PsgZgXHAtl33s771ChQldTCyLJ35dZrrT+fYOdfx2JbLIu+g4jeq8PBVx3Fzvo/n0jI7xjnvdN7FhtZ/wsiQfwQXcabXHs+hs2/AG1mp3JFjlC1YbCLijMAYk862sL71HwmstwvohtYrmVZzNLNqj+a0Rf+Pd7pW05F5FeERBJmoPmsSFNCYWEpz3XvxXa1exzhnT2Yjf9zxJaxPpTsjw87ux3i7/V/Zf9pnKiqTsIkYLFY2nBGoMO3dj7Kp/ft0pdZhSpM1CCzb7zgzo6XrfhY0foyYV8uCqS7nvmPis63rXoIBUuMYKTZ13MLCxk9X3pttEhsBV1SmgmzrvIP12z5BZ+oJAtox68YsRbFVKSMXeUs4HNVDLugmzJNWnEywk7fbvlk5gfJM4jgBZwQqRGBp3tp5NUbvH3ZvgMlIT3Fm1r2vApI5HJUjGZs+yF7DJ8c7HT8mk2sZ5Lgyk18TKGWrQpwRqBB7Mm/0qa4U4ilf2asngYmnWubWn0t9YlllhXQ4Rplcbhs+OfrWbCxMIGek2dJ+U0XlUhCUtFUjbk2gQsT9GQQDTO/EFeCbkcOjNr6EJTO/RkNyeYUldDhGH9+bSixKFpeNUhuqXwI5Y8euf2PmlFVMSVbC5bl6p3pKwRmBMmKWY1f3g3TsvodMduPeAnpGuPjrkSWgJ9lbIZ4MXzGWNX2PusRBlRXc4agQPgFS6JGTKBqhFeUWtTQtHTczpemHoy+UUTYjIGk/4KfAHMIJpJvN7HpJvwDy+f6nAW1mdqSkRcArwPpo3+NmdklZhCkRZwTKhFmat7edS3f6Oejj/gYQGMTxSFM86hfEfo2fdQbAUdUE1kmCLOmCQqeFhMYhNBTp3MYKCla2M2UJ0+o/K2kq8Iyk35nZufkDJH0baC/o85qZjVnNVWcEykTbrjvpTq+hmAGAqISdIEmAWc89lx8TeEoye2rRTBwOR9UwJXkMvupI2u5eqwL58XFhsfm4N9gicnkpV5xAlFRzS/S6U9IrhHVVXgZQ6Pt6DiOoBTxauIXhMmBm7Oy8ERjYpVOCWL6UhsCPNk/gezXMmfY5fK8yNVUdjrFias1J1CaWIcXCe18QE/iyAgMQfU+CvjWqRpFRcBGNpnreDTxR0HwCsNXM/ljQtljSc5JWSzphX1UZLm4kUAbaOq4hm3tryONiCiMTM3iA8DSVZPwAmhsuo7Hug6MvqMMxxkgeBzbfwfaOn7C1/TsEdIftBZOkHkYMw1PxKaOyYwa5kueDZkl6uuD9zQWVFvciqR74NfD3ZtZRsOt84OcF77cAC82sVdJy4G5Jh/bpM6o4I7CPBNZNx66b8QmiEJjBIx19gY/PtPoLaJ7+9UqI6HCMKzzVMLvx0wTZNbTvvqdgT+8nbd+rrZxQpT/lbzezFYMdIClOaABuN7N/L2iPAR8B9rr+WRgtmopePyPpNWApYcneiuCMwD6Sy71DmMItnFsLet3I/Q2CSJBMHEZT4xcrJKHDMU4JduIVjZ6J5qlzbZWTpXzeQQL+D/CKmX2nz+4PAOvMbGPB8U3ADjPLSToAOAh4vSzClIgzAvuI7zVD5M0QJxpZEj7TxPwFTJlyFuRSIBHzZ1GbXE5NYoWr9OWY9CTihxFPPYZZZu93RoQGwBPEKlVPIPThLtfZjgcuBF6Q9HzUdmVUV+U8ek8FAZxIWKkxS/jTcYmZ7SiXMKXgjMAAmKVIdd1KZvevQT6JuvNJ1J2H+sxTel4dU6dcTGfXjzHrjhaAQaqleeYPqEkOOnJ0OCYtDfUX09n1EyBT9IfIs+4KSWJg5fERNbNHGGBO2MwuKtL2a8KpozGjJCMg6Q2gk9BSZc1shaRrgdMJfSJfAz5mZv3Gb8X6lkf00cMsx67W88mlXyDv8dPdvoFs6g9MmdG/8MX0xivxvAY6dv0LQdBGPLaUGdO+6gyAwzEIsdg8Zk37Bi07/57C9QARjqq9Sg2WjeEsDFcdwxkJrDSz7QXvfwd8IaojfA3wBeCKEvuOCmaG5d4CDPn7j3jKJZv6A7nMy/R2+ewms2c12fQaYone1TYlj2kNlzGt4TLMzE31OBwlUlf3YWrb/4kg6EBSwfenlpraMyonyCROGzHiOAEze9BsbyL8x4EF5RFpZASZdaRa3k+q5VRSLatItZxEkF47onNlU4+DdRXZkyObfmrQvs4AOBylIyVpmHYdUi0QC78/qiMeP4S6KedXTpBJnEq61JGAAQ9KMuCmIn6xHwd+McK++4wFu0m1ngfWE4ltubdJ7biAmtmPIq9hWOeTPxuooV/wl+J43qx9F9jhcOyltu404omH6O66nVyuhWTNydTUrkIVK5VavT/wpVCqETjezDZLmg38TtI6M3sYQNIXCfNl3D7cvoVI+iTwSYCFCxcOS4ncnt9AsUUk20O27TJiDVeh2JKSz5eoPZM9nd/qV+tFxIjX/sWwZHM4HEMTiy1mauNVY/PhBlRpmuhSKGk6yMw2R3+3AXcBxwBI+hvgQ8AFZsVN6UB9ixx3s5mtMLMVTU1NJStgwQ6CzmsonrMnS5B6BNv+EYLdAw1U+uP5s6ifcSvyZoPqQLV4/kLqZ/0yGrY6HI6qwk0HDYykKYAXJUOaAnyQ0K/1VMKF4JPMbPdw+pZPfLDO6/CCLopXLSXKVbgHOv4Zq/kL5JXmexxLHkND81ME2fVADC+2xM33OxxVybDSRlQdpUwHNQN3RT+AMeBnZvZbSa8CScIpHojyYEuaB/yrma0aqG9ZNdjzACKISlP0ttSK/oVvYpD6L6hdVfKpJQ8/7qp7ORxVjYGVKU5gIjKkETCz14EjirQXnWSPpn9WDda3rCiOJOLmkyMgiILQPTx8vIKnd4ESoyqKw+GYoJQvYnjCMfFTSdeeDSSQREw+CcVJKE5Mfp/pG4Pke8dKSofDMZ6ZxGsCE94IqP4SSCwnjDHsSxw0BVSHpv0AqabS4jkcjvGOWegdVMpWhUz43EFSEs24BcusxdJPQ64FYgdC/FCUeRFUC8mVyJsy1qI6HI7xSpU+5ZfChDcCeRR/F4q/q3ejW9R1OBxDYlhuIP/C6qdqjIDD4XCMiPKmkp5wOCPgcDgczkXU4XA4JicGmBsJOBwOxyTFyldUZiLijIDD4Zj0TOaFYQ2Q921MkdQCvDnWcgzCLGDUi+SMIdWuHzgdq4WDzWzqvpxA0m8J/69KYbuZnbovnzfeGJdGYLwj6emJUCZzpFS7fuB0rBYmg46jzYSPGHY4HA7HyHFGwOFwOCYxzgiMjLKXyBxnVLt+4HSsFiaDjqOKWxNwOByOSYwbCTgcDsckxhmBQZB0tqSXJAWSVhS0XyDp+YItkHRktG+5pBckvSrpexrnNSkH0jHa9y5Jj0X7X1CUi7tadJS0SFJ3wXX8YcG+qtCxYP9CSbskXV7QVhU6Sjqm4BqukXRmwb4JpeOYYGZuG2ADlgEHA38AVgxwzOHA6wXvnwTeAwi4HzhtrPUYiY6EgYRrgSOi9zMBv8p0XAS8OECfqtCxYP+vgV8Cl1ebjkAdEItezwW2FbyfUDqOxeZGAoNgZq+Y2fohDjsf+DmApLlAg5k9ZuEd+FPgjFEWc58YRMcPAmvNbE10XKuZ5apMx6JUm46SzgBeB14qaKsaHc1st5llo7c1hOmAJqSOY4EzAvvOuURGAJgPbCzYtzFqm4gsBUzSA5KelfSPUXs16QiwWNJzklZLOiFqqxodJU0BrgC+0mdX1egIIOlYSS8BLwCXREahqnQcLSZ97iBJvwfmFNn1RTP7jyH6HgvsNrMX801FDhtz96sR6hgD3gscDewGHpL0DNBR5NiJquMWYKGZtUpaDtwt6VCq6zp+Bfiume3qMx1eTTpiZk8Ah0paBtwi6X7GqY7jjUlvBMzsA/vQ/Tx6RgEQPmksKHi/ANi8D+cvCyPUcSOw2sy2A0i6DzgKuI0q0dHMUkAqev2MpNcIR0DVdB2PBc6S9E1gGhBI2kO4RlAtOhb2f0VSF3AY4/Q6jjfcdNAIkeQBZwN35NvMbAvQKem4yAvho8Cgo4lxzAPAuyTVSYoBJwEvV5OOkpok+dHrA4CDCBf5q0ZHMzvBzBaZ2SLgOuDrZnZDNekoaXF0jyJpf8LF4zeqScdRZaxXpsfzBpxJ+DSRArYCDxTsex/weJE+K4AKqg2gAAAAh0lEQVQXgdeAG4gC8sbrNoSOf024mPgi8M1q0xH4y0i/NcCzwOnVpmOfY75Mb++gqtARuDC6js9H1/GMiarjWGwuYtjhcDgmMW46yOFwOCYxzgg4HA7HJMYZAYfD4ZjEOCPgcDgckxhnBBwOh2MS44yAw+FwTGKcEXA4HI5JjDMCDofDMYn5/ybOv0DZ2ITeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(glon,glat,c=area)\n",
    "cbar = plt.colorbar();\n",
    "cbar.set_label(\"area of each grid cell (km$^2$)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly, using a disk load with the same area isn't great if we are working with a lat/lon grid in high latitude regions because there are large variations in grid cell area accross the domain. Use a polar stereo grid for better results! Let's stick with this for now, though, and generate the 'G' matrix in our G.m=d formulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get locations of GPS stations \n",
    "gps_station_locations = sio.loadmat('data/gps_station_locations.mat')\n",
    "slon = gps_station_locations['slon']\n",
    "slat = gps_station_locations['slat']\n",
    "G = np.zeros((np.size(slat),np.size(glat)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geopy.distance import geodesic #Or use your favourite way to get distances between two lat/lons!\n",
    "fi = interp1d(dg_cm,ug_cm,kind ='slinear',fill_value=\"extrapolate\")\n",
    "for i in np.arange(0,len(slat)):\n",
    "    for j in np.arange(0,len(glon)):\n",
    "        dist = geodesic((slat[i],slon[i]),(glat[j],glon[j])).km\n",
    "        G[i,j] = (area[j]/area_disk)*fi(dist)\n",
    "sio.savemat('data/greens_matrix_025_deg_ALASKA.mat',{'slat':slat, 'slon':slon, 'glat':glat, 'glon':glon, 'G':G, 'grid_cell':grid_cell})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This 'G' matrix now can be used to estimate displacements due to a load (simply multiply G.m) or it can be used for inversions of loads using observed displacements (solve the matrix equation G.m = d)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
